快速傅里叶变换模块(fft)
法国科学家傅里叶提出，任何一条周期曲线，无论多么跳跃或不规则，都能表示成一组光滑正弦曲线叠加之和。
傅里叶变换的目的是可将时域（即时间域）上的信号转变为频域（即频率域）上的信号，随着域的不同，对同一个事物的了解角度也就随之改变，因此在时域中某些不好处理的地方，在频域就可以较为简单的处理。这就可以大量减少处理信号存储量。
例如：弹钢琴
假设有一时间域函数：y = f(x)，根据傅里叶的理论它可以被分解为一系列正弦函数的叠加，他们的振幅A，频率ω或初相位φ不同：
所以傅里叶变换可以把一个比较复杂的函数转换为多个简单函数的叠加，看问题的角度也从时间域转到了频率域，有些的问题处理起来就会比较简单。
导入快速傅里叶变换所需模块
import numpy.fft as nf
freqs = np.fft.fftfreq(采样数量, 采样周期)
np.fft.fft(原函数值序列) -> 目标函数值序列(复数)
np.fft.ifft(目标函数值序列(复数))->原函数值序列
案例：针对合成波做快速傅里叶变换，得到一组复数序列；再针对该复数序列做逆向傅里叶变换得到新的合成波并绘制。
ffts = nf.fft(sigs6) sigs7 = nf.ifft(ffts).real mp.plot(times, sigs7, label=r'$\omega$='+str(round(1 / (2 * np.pi),3)), alpha=0.5, linewidth=6)
案例：针对合成波做快速傅里叶变换，得到分解波数组的频率、振幅、初相位数组，并绘制频域图像。
# 得到分解波的频率序列 freqs = nf.fftfreq(times.size, times[1] - times[0]) # 复数的模为信号的振幅（能量大小） ffts = nf.fft(sigs6) pows = np.abs(ffts) mp.subplot(122) mp.title('Frequency Domain', fontsize=16) mp.xlabel('Frequency', fontsize=12) mp.ylabel('Power', fontsize=12) mp.tick_params(labelsize=10) mp.grid(linestyle=':') mp.plot(freqs[freqs >= 0], pows[freqs >= 0], c='orangered', label='Frequency Spectrum') mp.legend() mp.tight_layout() mp.show()
""" 傅里叶变换拆方波 """ import numpy as np import matplotlib.pyplot as mp import numpy.fft as nf x = np.linspace(-2*np.pi, 2*np.pi, 1000) y = np.zeros(x.size) for i in range(1, 1000): y += 4*np.pi/(2*i-1) * np.sin((2*i-1)*x) mp.figure('FFT',facecolor='lightgray') mp.subplot(121) mp.title('Time Domain',fontsize=18) mp.grid(linestyle=":") mp.plot(x, y, label=r'$y$') #针对方波y,做fft comp_ary = nf.fft(y) y2 = nf.ifft(comp_ary) mp.plot(x,y2,color = 'orangered',linewidth=5,alpha=0.5,label='$y$') #绘制频域图像(频率/能量) mp.subplot(122) freqs = nf.fftfreq(y.size,x[1]-x[0]) pows = np.abs(comp_ary)#负数的模 mp.title('Frequency Domain',fontsize=18) mp.grid(linestyle=":") mp.plot(freqs[freqs>0],pows[freqs>0],color= 'orangered',label='frequency') mp.legend() mp.show()
基于傅里叶变换的频域滤波
含噪信号是高能信号与低能噪声叠加的信号，可以通过傅里叶变换的频域滤波实现降噪。
通过FFT使含噪信号转换为含噪频谱，去除低能噪声，留下高能频谱后再通过IFFT留下高能信号。
案例：基于傅里叶变换的频域滤波为音频文件去除噪声。
读取音频文件，获取音频文件基本信息：采样个数，采样周期，与每个采样的声音信号值。绘制音频时域的：时间/位移图像。
import numpy as np import numpy.fft as nf import scipy.io.wavfile as wf import matplotlib.pyplot as mp sample_rate, noised_sigs = wf.read('../data/noised.wav') noised_sigs = noised_sigs / 2 ** 15 times = np.arange(len(noised_sigs)) / sample_rate mp.figure('Filter', facecolor='lightgray') mp.subplot(221) mp.title('Time Domain', fontsize=16) mp.ylabel('Signal', fontsize=12) mp.tick_params(labelsize=10) mp.grid(linestyle=':') mp.plot(times[:178], noised_sigs[:178],c='orangered', label='Noised') mp.legend() mp.show()
基于傅里叶变换，获取音频频域信息，绘制音频频域的：频率/能量图像。
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=10) mp.grid(linestyle=':') mp.semilogy(freqs[freqs >= 0],noised_pows[freqs >= 0], c='limegreen',label='Noised') mp.legend()
将低能噪声去除后绘制音频频域的：频率/能量图像。
fund_freq = freqs[noised_pows.argmax()] noised_indices = np.where(freqs != fund_freq) filter_ffts = noised_ffts.copy() filter_ffts[noised_indices] = 0 filter_pows = np.abs(filter_ffts) mp.subplot(224) mp.xlabel('Frequency', fontsize=12) mp.ylabel('Power', fontsize=12) mp.tick_params(labelsize=10) mp.grid(linestyle=':') mp.plot(freqs[freqs >= 0], filter_pows[freqs >= 0],c='dodgerblue', label='Filter') mp.legend()
基于逆向傅里叶变换，生成新的音频信号，绘制音频时域的：时间/位移图像。
filter_sigs = nf.ifft(filter_ffts).real mp.subplot(223) mp.xlabel('Time', fontsize=12) mp.ylabel('Signal', fontsize=12) mp.tick_params(labelsize=10) mp.grid(linestyle=':') mp.plot(times[:178], filter_sigs[:178],c='hotpink', label='Filter') mp.legend()
重新生成音频文件。
wf.write('../../data/filter.wav',sample_rate,(filter_sigs * 2 ** 15).astype(np.int16))
""" 频域滤波 """ import numpy as np import numpy.fft as nf import matplotlib.pyplot as mp import scipy.io.wavfile as wf # 采样率,采样位移 sample_rate, sigs = wf.read('noised.wav') print(sample_rate) # 44100 print(sigs.shape) # (220500,) sigs = sigs/(2**15) times = np.arange(sigs.size) / sample_rate print(times) """ [0.00000000e+00 2.26757370e-05 4.53514739e-05 ... 4.99993197e+00 4.99995465e+00 4.99997732e+00] """ # 绘图 mp.figure('Filter', facecolor='lightgray') mp.subplot(221) mp.title('Time Domain', fontsize=16) mp.ylabel('Signal', fontsize=14) mp.grid(linestyle=':') mp.plot(times[:178], sigs[:178], color='dodgerblue', label='Noised Signal') mp.legend() # fft变换后绘制频域图(频率/能量) comp_ary = nf.fft(sigs) pows = np.abs(comp_ary) freqs = nf.fftfreq(sigs.size, times[1] - times[0]) mp.subplot(222) mp.title('Frequency Domain', fontsize=16) mp.ylabel('pow', fontsize=14) mp.grid(linestyle=":") mp.semilogy(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='Freq') mp.legend() # 在频谱中去除噪声 maxpow_freq = freqs[pows.argmax()] noised_inds = np.where(freqs != maxpow_freq) # 噪声的复数全抹为0 comp_ary[noised_inds] = 0 # 绘制降噪之后的频谱图像 pows = np.abs(comp_ary) mp.subplot(224) mp.ylabel('pow', fontsize=14) mp.grid(linestyle=':') mp.plot(freqs[freqs > 0], pows[freqs > 0], color='orangered', label='Freq') mp.legend() # 逆向傅里叶变换,输出时域函数图像 filter_sigs = nf.ifft(comp_ary) mp.subplot(223) mp.ylabel('Signal', fontsize=14) mp.grid(linestyle=":") mp.plot(times[:178], filter_sigs[:178], color='dodgerblue', label='Filter Signal') mp.legend() # 保存文件 wf.write('filter.wav', sample_rate, (filter_sigs * 2 ** 15).astype(np.int16)) mp.tight_layout() mp.show()
