文章目录
- 1. 傅里叶变换
- 2. 离散傅里叶变换
- 3. 快速傅里叶变换推导过程
说起傅里叶变换和逆变换,我的内心是抗拒的,但其实没必要深究其原理,会用就行
在计算机的世界,只能处理离散的信号,数字图像处理也是处理离散的图像数据.这里会重点介绍下FFT(快速傅里叶变换).
1. 傅里叶变换
- 傅里叶变换的正变换和逆变换公式
傅里叶变换是我们既爱又恨的公式,"恨"的是,它不太容易理解.爱的是它异常强大.关于它的推导我就不推了,关键是我也不会,公式记住即可.
下面是用AI生成的一个傅里叶变换的例子,用来将3种频率分别为5,15,30的周期信号,先进行叠加,然后再用傅里叶变换解调出各中频率的例子.
由于计算机只能处理离散的信号,下面再脚本一开始的地方,可以看到先对图像进行采样,然后做快速傅里叶变换.
importnumpyasnpimportmatplotlib.pyplotasplt# --- 1. 定义参数 ---# 采样参数sampling_rate=1000# 采样频率 (Hz)dt=1/sampling_rate# 时间间隔 (s)duration=2.0# 信号持续时间 (s)N=int(duration*sampling_rate)# 采样点总数# 时间轴t=np.linspace(0,duration,N,endpoint=False)# 输入信号的频率成分 (Hz)f1,f2,f3=5,15,30# --- 2. 生成各种输入信号 ---# 正弦波信号 s1(t) = sin(2π * f1 * t)s1=np.sin(2*np.pi*f1*t)# 余弦波信号 s2(t) = 0.5 * cos(2π * f2 * t)s2=0.5*np.cos(2*np.pi*f2*t)# 正弦波信号 s3(t) = 0.3 * sin(2π * f3 * t)s3=0.3*np.sin(2*np.pi*f3*t)# --- 3. 生成叠加信号 ---# 叠加信号 x(t) = s1(t) + s2(t) + s3(t)combined_signal=s1+s2+s3# --- 4. 对叠加信号进行快速傅里叶变换 (FFT) ---# 计算FFTfft_result=np.fft.fft(combined_signal)# 计算对应的频率轴# np.fft.fftfreq(n, d=dt) 计算FFT后每个点对应的频率# n 是FFT的点数，d是采样时间间隔freqs=np.fft.fftfreq(N,d=dt)# --- 5. 为了绘图，只取频谱的前半部分（正频率部分） ---# FFT结果是对称的，通常只关心正频率部分# 幅度谱取绝对值magnitude_spectrum=np.abs(fft_result)half_N=N//2positive_freqs=freqs[:half_N]positive_magnitude_spectrum=magnitude_spectrum[:half_N]# --- 6. 绘制图形 ---# 创建一个包含3个子图的图形，采用垂直堆叠（1列，3行）# figsize 控制整体图形大小plt.figure(figsize=(12,10))# 第一个子图：绘制原始的几个信号plt.subplot(3,1,1)# (总行数, 总列数, 当前子图序号)plt.plot(t,s1,label=f'信号 1:{f1}Hz',alpha=0.7)plt.plot(t,s2,label=f'信号 2:{f2}Hz',alpha=0.7)plt.plot(t,s3,label=f'信号 3:{f3}Hz',alpha=0.7)plt.title('原始信号')plt.xlabel('时间 (s)')plt.ylabel('振幅')plt.legend()plt.grid(True,linestyle='--',alpha=0.6)# 第二个子图：绘制叠加后的信号plt.subplot(3,1,2)plt.plot(t,combined_signal,label='叠加信号',color='purple')plt.title('叠加信号')plt.xlabel('时间 (s)')plt.ylabel('振幅')plt.legend()plt.grid(True,linestyle='--',alpha=0.6)# 第三个子图：绘制FFT结果的频域图plt.subplot(3,1,3)plt.plot(positive_freqs,positive_magnitude_spectrum,label='频谱',color='red')plt.title('傅里叶变换 (频域)')plt.xlabel('频率 (Hz)')plt.ylabel('幅度')plt.xlim(0,50)# 限制x轴范围，只看低频部分，更容易观察到目标频率plt.grid(True,linestyle='--',alpha=0.6)# 在频域图中标注出峰值，即原始频率# 找到幅度最大的几个点的索引peaks_indices=np.argsort(positive_magnitude_spectrum)[-10:]# 取最大的10个点的索引# 获取这些点对应的频率值和幅度值peak_freqs=positive_freqs[peaks_indices]peak_magnitudes=positive_magnitude_spectrum[peaks_indices]# 只标注接近我们设定频率的峰值foriinrange(len(peak_freqs)):ifabs(peak_freqs[i]-f1)<0.5orabs(peak_freqs[i]-f2)<0.5orabs(peak_freqs[i]-f3)<0.5:plt.annotate(f'{peak_freqs[i]:.1f}Hz',xy=(peak_freqs[i],peak_magnitudes[i]),xytext=(peak_freqs[i]+2,peak_magnitudes[i]+max(positive_magnitude_spectrum)*0.02),arrowprops=dict(arrowstyle='->',color='black'))plt.tight_layout()# 自动调整子图间距，防止重叠plt.show()# 显示图形- 实验结果图
下图是上下结构,最上面的原始信号,中间的是叠加后的信号,最后是用离散傅里叶变换获取到的频谱图,可以看到信号的频率和一开始的输入信号是一致的.
2. 离散傅里叶变换
记住这个公式
3. 快速傅里叶变换推导过程
离散傅里叶变换已成为数字信号处理的重要工具,但是它的计算量很大,计算时间长,在实时图像的分析上显然性能不够.所以为了解决这个矛盾 美国的两位科学家 詹姆斯·库利（James Cooley） 和 约翰·图基（John Tukey）提出了快速傅里叶变换的方法
快速傅里叶变换以数字序列长度N的状况分为以下3种算法
- N为2的整数幂的算法
- N为高复合数的算法
- N为素数的算法(素数:大于1的自然数中，除了1和它本身以外不再有其他因数的自然数)
下面是我自己以N=8的数字序列,参考教材自己也推导了一遍.这里记录下
从上面的推导过程可以知道FFT就是一个分而治之,递归的方法.上面的4点FFT前面,实际运算时还有个2点FFT,图中没有列出.
如果是32点的数字序列,该如何计算离散傅里叶变换.利用上面的方法,它最后一定是下面这种小的FFT叠加在一起的. 如下图,这只是一个示意图,一些加权系数没有标在图中,大概知道什么意思了,主打一个会用
2 4 8 16
我的理解是,按着分而治之的方法,逐步将复杂的FFT分解成简单简单的数值运算.如下所示
