一,前期回顾:
上一节主要讲了傅里叶基数到傅里叶变换,其主要的思想可以总结为两句话。
- 对于傅里叶基数,“一个周期连续的波形可以由多个其周期整数倍的波形组合而成!” 由此给出了公式并进行了系数推导。
- 从傅里叶基数到傅里叶变换的 idea 为 “ 把周期信号的周期逐渐扩大,当接近无穷大,这样周期信号不就成了非周期的!”
最终我们得出的结论是:我们能够将时域中非周期连续的信号通过傅里叶变换华丽转身到频域中连续非周期的频谱。这在理论上实现了巨大的成功,但是又回到实际中,我们的PC 不能处理连续的信号,数字信号才是PC的菜啊!
本文主要围绕着怎样对信号进行处理,可以使PC能handle!
二,傅里叶变换到离散时间傅里叶变换 (FT 到 DTFT)
首先,我们要明确一点就是傅里叶变换处理的是连续的系统,连续系统,连续系统,重要事情说三遍!
很直接想到的一个idea就是,采样啊,让狄拉克梳状函数上啊!
我们也有一篇专门讲解采样的文章!
马尚先生:采样定理,频谱混叠和傅里叶变换 深入理解zhuanlan.zhihu.com
千万别走偏,精彩继续,我们继续往下看。
冲击采样序列即狄拉克梳状函数的表达式为:
采样后的信号为:
形象点如图:
Fig.1 狄拉克梳状函数对连续信号进行采样过程简图
连续信号的傅里叶变换为:
采样后信号的傅里叶变换为:
改变下求和积分的顺序,并利用冲击函数的筛选性质。
。可以得到离散时间傅里叶变换公式为:
进一步简化对于数字信号而言,我们只有x(t) 采样后的信号x[n],第n个采样发生在t=nTs 时因此可以将连续信号的傅里叶变换改写为:
小结:
到了这里我们在来总结下上面所做的工作。此时我们已经完成了离散时间傅里叶变换DTFT的公式的推导。在傅里叶变换的基础上,我们对时域的信号进行了采样,将时域连续的信号变为了离散信号并进行傅里叶变换,我们神奇的发现在频域中还是连续的谱线。
因为时域中采样,相当于在频域中以采样频率对信号频率进行周期延拓。
Fig.2 解释时域和频域对偶问题
问题:
此时时域信号离散了,但是频域信号还是连续的,我们要的是时域和频域都是离散的这样PC才可能处理,现在又该怎么办呐?
下面是我自己的理解,思路是这样的!
现在只有频域是连续的,可以采用和时域相同的方法,再对频域进行采样不就可以了吗!在利用 “一个域与信号相乘等同于在另一个域中卷积”原理。我们在频域中将连续的频谱与狄拉克梳状函数相乘,这样时域等同于把信号以采样频率为周期进行延拓。
这样做的前提是提前在时域中选取N个点作为一个周期,这样在延拓的时候才有东西可以延拓。
思路是这样下面我们来具体的实现下!
三,离散傅里叶变换 DFT
首先在时域中选取N个点,说白了就是对时域连续信号x(t) 进行N点采样,然后将N点采样信号进行周期延拓,虚拟成周期离散的信号并将其进行离散傅里叶变换,得到的频域谱图也是周期离散的。有的小同志就发现了你时域中信号周期延拓,频域也是周期的,是不是只观察频域中的一个区间的取值就可以知道信号的全貌。对的,一般我们只取主值区间[0,N-1]即可,因为其他都是相同的。
有了思路我们在理论上推导下。
对连续的信号进行N点的采样得到信号的表达式为,这个公式似曾相识啊! 注意上下标,哈哈
将假想为N个点做周期延拓,取主值区间进行傅里叶变换,有:
积分和求和的顺序调换下,并利用狄卡拉梳状函数的筛选特性得:
当把T换为N时,即只有在 t=nTs 时候有值,上式可以化简为
此式为离散周期傅里叶变换,理论上讲周期信号的频谱为无穷多的,但是由于
的周期性,即
一般只取0 到N-1. 此式即为离散傅里叶公式。
有的同志又注意到了,最后的n 是对应的点数而不是频率啊,这怎样转换到频率那?看下面的公式:
也可以写为
, 其中N为做FFT的点数,fs 为采样频率。
所以离散周期傅里叶变换后的第k个数对应的频率是
,fs/N 可以理解为频率分辨率。
四,实际应用中的小问题
1,为什么最好进行2次幂个点数的FFT运算?
在实际的应用中如果直接按DFT变换进行计算,当序列长度N很大时,计算量会非常大,所需时间也很长,因此常用的是DFT的一种快速计算算法快速傅里叶变换FFT。
最常用的FFT算法是基于时间抽取的基2-FFT算法和基于频率抽取的基2-FFT算法,这种算法的特点在于FFT会把一次大的DFT分割成几个小的DFT,这样递归式地细分下去,例如有8个采样点的FFT,首先会把最外层的8点运算分成两个4点FFT的奇偶组合,第二层FFT又分成四个两点FFT的奇偶组合,并且由此计算出的频谱中很有趣的一点在于对于实数输出的数组,后面一半和前面一半正好对称相同,对于虚数输出的数组,后面一半是前面数组对称后乘上负1,因此,我们只需要算出FFT的一半即可求出全部。
所以在进行FFT的时候一般选取2 的次幂个点数进行快速傅里叶变换。
2,显示分辨率和实际分辨率。
在选取N个点数做FFT以后一遍要乘以频率分辨率来得到相应的频率左边,频率分辨率的公式为Fre = Fs/N, 这里N的取值不同,得到结果的意义也就不同,当除以信号本身的长度时候为实际的频率分辨率,当信号不够2 的整数次幂,进行补零凑齐以后,此时N如果取补零以后的个数,那么此时的分辨率为显示分辨率,实际分辨率并没有提高。
例如,要求的频率分辨率要达到0.5 Hz 此时的采样率 fs =10 kHz, 那么要满足要求,采样的点数必须最少为20 K 个点,否者无法满足要求。
3,MATLAB 中FFT 的实现
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
4,MATLAB 中的FFTshift 是怎么个意思?
信号以采样频率fs进行采样会在在 fs/2处混叠,所以实信号fft的结果中前半部分对应[0, fs/2],后半部分对应[ -fs/2, 0],大于fs/2的部分的频谱实际上是信号的负频率加fs的结果。故要得到正确的负频率的结果,只需将视在频率减去采样频率fs即可得到频谱对应的真实负频率。
进行N 点的FFT 最后输出的结果是以fs/2 为中心的频谱,如果只想显示正的频率,只显示[0, fs/2] 即可,如果想要显示负的频率,此时就需要fftshift 一下了,把 对称的中心从fs/2 移动到 0 点。具体的结果如下:
Fig MATLAB 中 fftshift 作用和实现
完整MATLAB 代码
clear
fs=200;N=512; %采样频率和数据点数
n=0:N-1;t=n/fs; %时间序列
x=0.9*sin(2*pi*20*t)+2*sin(2*pi*80*t); %信号
y1=fft(x,N); %对信号进行FFT
y2=fftshift(y1);
mag1=2*abs(y1)/N; %求得Fourier变换后的振幅
mag2=2*abs(y2)/N;
f1=n*fs/N; %频率序列
f2=n*fs/N-fs/2;
fig1=figure();
subplot(3,1,1);
plot(f1,mag1,'r'); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');
title('图1 正常FFT ');
grid on;
subplot(3,1,2),plot(f2,mag1,'b'); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');
title('图2 直接平移fs/2,频率不是很对');
grid on;
subplot(3,1,3);
plot(f2,mag2,'g'); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('图3 把频率对称点移到0点');
grid on;
saveas (fig1,'fftshift11.jpg')
python 中 FFT 实现方法
详见以下链接:
马尚先生:Python 中 FFT 快速傅里叶分析zhuanlan.zhihu.com