numpy的快速傅里叶变换

快速傅里叶变换(fft)

什么是傅里叶变换?

傅里叶定理: 任何一条周期曲线, 无论多么跳跃或不规则, 都能表示成一组光滑正弦曲线叠加之和. 傅里叶变换即是把这条周期曲线拆解成一组光滑正弦曲线的过程.

傅里叶变换的目的是将时域(时间域)上的信号转变为频域(频率域)上的信号, 随着域的不同,对同一个事物的了解角度也随之改变. 因此在时域中某些不好处理的地方, 放在频域中就可以较为简单的处理. 这样可以大量减少处理的数据量.

傅里叶定理:
$$y=A_{1}sin({\omega }_{1}x+{\varphi }_1)+A_{2}sin({\omega }_{2}x+{\varphi }_2)+..+C$$
快速傅里叶变换相关API:

import numpy.fft as nf
# 通过采样数与采样周期,得到fft分解所得曲线的频率数组
# 采样周期: x轴相邻两点的距离
freqs = nf.fftfreq(采样数量, 采样周期)

# 通过原函数值序列, 经过fft后, 得到复数数组
# 复数数组的长度即是拆解出正弦函数的个数
# 复数数组每个元素的模,代表每个正弦曲线的振幅
# 复数数组每个元素的辅角,代表每个正弦曲线的相位角
复数序列 = nf.fft(原函数值序列)

# 逆向傅里叶变换
原函数值序列 = nf.ifft(复数序列)

案例: 基于傅里叶变换, 拆解方波.

"""
demo04_fft.py  傅里叶变换
"""
import numpy as np
import matplotlib.pyplot as mp
import numpy.fft as nf

x = np.linspace(-2*np.pi, 2*np.pi, 1000)

# 叠加1000条曲线
y = np.zeros(x.size)
for i in range(1, 1000):
	y += 4/(2*i-1)*np.pi * np.sin((2*i-1)*x)

# 对y做傅里叶变换, 绘制频域图像
ffts = nf.fft(y)
# 获取傅里叶变换的频率序列
freqs = nf.fftfreq(x.size, x[1]-x[0])
pows = np.abs(ffts)

mp.figure('FFT', facecolor='lightgray')
mp.subplot(121)
mp.grid(linestyle=':')
mp.plot(x, y,linewidth=2)
mp.subplot(122)
mp.grid(linestyle=':')
mp.plot(freqs[freqs>0], pows[freqs>0], 
	    c='orangered')
# 通过复数数组,经过ifft操作, 得到原函数
y2 = nf.ifft(ffts)
mp.subplot(121)
mp.plot(x, y2, linewidth=7, alpha=0.5)

mp.show()

基于傅里叶变换的频域滤波

含噪信号是高能信号与低能噪声叠加的信号, 可以通过傅里叶变换的频域滤波实现简单降噪.

通过FFT使含噪信号转换为含噪频谱, 手动取出低能噪声, 留下高能频谱后,再通过IFFT生成高能信号.

1. 读取音频文件, 获取音频的基本信息: 采样个数/采样周期/每个采样点的声音值. 绘制音频的时域: 时间/位移图像.

sample_rate, noised_sigs = \
	wf.read('../da_data/noised.wav')
print(noised_sigs.shape) 
times = np.arange(len(noised_sigs))/sample_rate

mp.figure('Filter', facecolor='lightgray')
mp.subplot(2,2,1)
mp.title('Time Domain', fontsize=16)
mp.ylabel('Signal', fontsize=12)
mp.tick_params(labelsize=8)
mp.grid(linestyle=':')
mp.plot(times[:178], noised_sigs[:178], 
	c='dodgerblue', label='Noised Sigs')
mp.legend()
mp.show()

2. 基于傅里叶变换, 获取音频频域信息, 绘制频域: 频率/能量图像.

freqs = nf.fftfreq(times.size, 1/sample_rate)
noised_ffts = nf.fft(noised_sigs)
noised_pows = np.abs(noised_ffts)
mp.subplot(222)
mp.title('Frequency Domain', fontsize=16)
mp.ylabel('Power', fontsize=12)
mp.tick_params(labelsize=8)
mp.grid(linestyle=":")
mp.semilogy(freqs[freqs>0], noised_pows[freqs>0],
	c='orangered', label='Noised')
mp.legend()

3. 将低频噪声去除后绘制音频频域: 频率/能量图像.

noised_inds = np.where(freqs != fund_freq)
filter_ffts = noised_ffts.copy()
filter_ffts[noised_inds] = 0
filter_pows = np.abs(filter_ffts)
mp.subplot(224)
mp.title('Frequency Domain', fontsize=16)
mp.ylabel('Power', fontsize=12)
mp.tick_params(labelsize=8)
mp.grid(linestyle=":")
mp.plot(freqs[freqs>0], filter_pows[freqs>0],
	c='orangered', label='Filter')
mp.legend()

4. 基于逆向傅里叶变换,生成时域的音频信号, 绘制时域: 时间/位移图像.

filter_sigs = nf.ifft(filter_ffts)
mp.subplot(2,2,3)
mp.title('Time Domain', fontsize=16)
mp.ylabel('Signal', fontsize=12)
mp.tick_params(labelsize=8)
mp.grid(linestyle=':')
mp.plot(times[:178], filter_sigs[:178], 
	c='dodgerblue', label='Filter Sigs')
mp.legend()

5. 生成音频文件

wf.write('../da_data/filter.wav', sample_rate, filter_sigs.astype(np.int16))