背景
晚上看到一篇文章，介绍了【对20世纪科学与工程的发展和实践影响最大的10大算法】，看到了 FFT 算法，于是写点。
Fast Fourier Transform (FFT) 是一种用于计算离散傅里叶变换 (DFT) 的算法。
傅里叶变换是一种将信号从时间域转换到频率域的数学变换，而FFT通过优化计算过程使这一变换更加高效。
快速傅里叶变换
考虑以下场景，如何高效的解决？
Q：假设你有一段音乐，它由多个频率成分组成。
A：FFT可以帮助你找到这些成分的频率和强度。
比如说，你可以通过FFT发现这段音乐中主要的频率成分是100Hz、200Hz和300Hz，且它们的强度分别是50、30和20。
什么是傅里叶变换？
傅里叶变换将复杂的信号分解成一系列简单的正弦波（即不同频率的正弦和余弦波）。
这类似于将一首复杂的交响乐分解成单独的乐器声部，使我们能够分析每种频率成分的强度。
什么是快速傅里叶变换？
FFT，英文全称，Fast Fourier Transform。
傅里叶变换的直接计算（即DFT）非常耗时，特别是当信号长度较长时。
FFT是一种优化算法，使傅里叶变换的计算变得更快。
它利用信号的对称性和周期性，将计算复杂度从 O ( N 2) 降低到 O ( N log N ) 。
电力系统的应用场景
一个完整的可视化示例
这是一个使用快速傅立叶变换（FFT）进行音频噪声消除的示例，用于模拟和实现音频信号的降噪系统。
1、导入包
import matplotlib.pyplot as pltimport numpy as np
2、创建一个包含多个频率成分的信号，并为信号添加了噪声
# 设置图形大小和字体大小plt.rcParams['figure.figsize']=[16,12]plt.rcParams.update({'font.size':18})# 创建时间轴dt =0.001t = np.arange(0,1, dt)# 创建包含10种不同频率的信号f1 = np.sin(2* np.pi *20* t)f2 = np.sin(2* np.pi *50* t)f3 = np.sin(2* np.pi *70* t)f4 = np.sin(2* np.pi *90* t)f5 = np.sin(2* np.pi *120* t)f6 = np.sin(2* np.pi *150* t)f7 = np.sin(2* np.pi *180* t)f8 = np.sin(2* np.pi *210* t)f9 = np.sin(2* np.pi *240* t)f10 = np.sin(2* np.pi *270* t)# 合并多个频率f_clean = f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 + f10# 为信号添加噪声f_noisy = f_clean +4* np.random.randn(len(t))
3、生成干净信号的音频形式
# 导入音频播放功能from IPython.display import Audio# 扩展信号f_clean_extended = np.tile(f_clean, 100)# 转换信号数据类型f_clean_audio_data = f_clean_extended.astype(np.float32)# 播放音频display(Audio(f_clean_audio_data , autoplay=True, rate=44100))
4、将包含白噪声的信号转换为音频形式
# 将信号转换为音频数据f_noisy_extended = np.tile(f_noisy, 100)f_noisy_audio_data = f_noisy_extended.astype(np.float32)display(Audio(f_noisy_audio_data, autoplay=True, rate=44100))
5、绘制干净信号和包含白噪声的信号的时域波形图
# 绘制干净和噪声信号的时域波形图plt.plot(t , f_noisy, color='c', linewidth=1.5, label='Noisy')plt.plot(t, f_clean, color='k', linewidth=2, label='Clean')# 设置绘图限制和标签plt.xlim(t[0], t[-1]) # 设置x轴的范围为从t的第一个值到最后一个值plt.legend() #显示图例plt.show()
6、进行频谱分析，并绘制包含噪声信号的时域和频域图
n = len(t)fhat = np.fft.fft(f_noisy,n)PSD = fhat * np.conj(fhat)/ n #power spectral densityfreq =(1/(dt*n))* np.arange(n)L = np.arange(1,np.floor(n/2), dtype ='int')fig,axs = plt.subplots (2,1)plt.sca(axs[0])plt.plot(t,f_noisy,color='c', label='Noisy')plt.plot(t,f_clean,color='k', label='Clean')plt.xlim(t[0],t[-1])plt.legend()plt.sca(axs[1])plt.plot(freq[L],PSD[L],color='c', label='Noisy')plt.xlim(freq[L[0]], freq[L[-1]])plt.legend()plt.show()
7、利用功率谱密度（PSD）清除信号中的噪声成分
# 使用功率谱密度（PSD）清除噪声indices = PSD > 150 # 找出所有功率大的频率PSDclean = PSD * indices # 将其他频率清零fhat = indices * fhat # 将Y中小的Fourier系数清零ffilt = np.fft.ifft(fhat) # 反FFT
# 将信号转换为音频数据ffilt_extended = np.tile(ffilt , 100)ffilt_audio_data = ffilt_extended.astype(np.float32)display(Audio(ffilt_audio_data, autoplay=True, rate=44100))
8、绘制多个子图，展示不同信号处理步骤的结果
# 绘制多个子图fig,axs = plt.subplots(3,1)plt.sca(axs[0])plt.plot(t, f_noisy, color='c', label='Noisy')# 绘制包含噪声的信号时域图plt.plot(t, f_clean, color='k', label='Clean')# 绘制干净信号时域图plt.xlim(t[0], t[-1])# 设置x轴的范围为从时间数组t的第一个值到最后一个值plt.legend()# 显示图例plt.sca(axs[1])plt.plot(t, ffilt, color='k', label='Filtered')# 绘制滤波后的信号时域图plt.xlim(t[0], t[-1])# 设置x轴的范围为从时间数组t的第一个值到最后一个值plt.legend()# 显示图例plt.sca(axs[2])plt.plot(freq[L], PSD[L], color='c', label='Noisy')# 绘制包含噪声信号的频域图plt.plot(freq[L],PSDclean[L], color='k', label='Filtered')# 绘制经过滤波后信号的频域图plt.xlim(freq[L[0]], freq[L[-1]])# 设置x轴的范围为从频率数组L的第一个值到最后一个值plt.legend()# 显示图例plt.show() # 显示绘制的图形
也可以参看之前文章： 傅里叶变换及Python实现
