加载中…MATLAB中使用FFT做频谱分析时频率分辨率问题
(2010-12-25 22:23:44)
标签：
杂谈 |
分类： FFT |
MATLAB中使用FFT做频谱分析时频率分辨率问题
最近做FFT时，使用的采样频率和信号长度的取舍一直没有搞清楚，后来在论坛上发了一个贴子《总结一下使用FFT和维纳-辛钦定理求解PSD问题》（讨论见http://www.ilovematlab.cn/thread-27150-1-1.html，特别感谢会员songzy41，他的问题给了我很大启示），跟帖中给了我不少启示，并且让我对“频率分辨率”这个概念有了更深入的理解。再次一并感谢论坛的高手们。
频率分辨率，顾名思义，就是将信号中两个靠的很近的频谱分开的能力。
信号x(t)长度为Ts，通过傅氏变换后得到X，其频率分辨率为Δf=1/T（Hz），若经过采样后，假设采样频率为fs=1/Ts，而进行频谱分析时要将这个无穷长的序列使用窗函数截断处理，假设使用矩形窗，我们知道，矩形窗的频谱为sinc函数，主瓣宽度可以定义为2*pi/M，M为窗宽，那么，时域相乘相当于频域卷积，频域内，这一窗函数能够分辨出的最近频率肯定不可能小于2*pi/M了，也就是如果数据长度不能满足2*pi/M<|w2-w1|（w2,w1为两个靠的很近的频率），那么在频谱分析时，频谱上将不能分辨出这两个谱，由于w2-w1=2*pi(f2-f1)/fs=2*pi*Δf/fs也就是2*pi/M<2*piΔf/fs，得到Δf的限制为fs/M，这就是窗函数宽度的最小选择，就是说，根据Shannon采样定理确定了采样频率后，要根据靠的最近的谱峰来确定最小的采样长度，这样，所作出来的频谱才能分辨出那两个谱峰，也就是拥有了相应的频率分辨率。
几个例子：
[CODE]
Fs = 80;
n = 0:1/Fs:1023*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
figure(1);
subplot(211)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长1024，包含整周期');
n = 0:1/Fs:979*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
subplot(212)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长980，包含整周期');
n = 0:1/Fs:399*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
figure(2);
subplot(211)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长400，不包含整周期');
Fs = 20;
n = 0:1/Fs:1024*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
subplot(212)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长1024，采样率过小');
最近做FFT时，使用的采样频率和信号长度的取舍一直没有搞清楚，后来在论坛上发了一个贴子《总结一下使用FFT和维纳-辛钦定理求解PSD问题》（讨论见http://www.ilovematlab.cn/thread-27150-1-1.html，特别感谢会员songzy41，他的问题给了我很大启示），跟帖中给了我不少启示，并且让我对“频率分辨率”这个概念有了更深入的理解。再次一并感谢论坛的高手们。
频率分辨率，顾名思义，就是将信号中两个靠的很近的频谱分开的能力。
信号x(t)长度为Ts，通过傅氏变换后得到X，其频率分辨率为Δf=1/T（Hz），若经过采样后，假设采样频率为fs=1/Ts，而进行频谱分析时要将这个无穷长的序列使用窗函数截断处理，假设使用矩形窗，我们知道，矩形窗的频谱为sinc函数，主瓣宽度可以定义为2*pi/M，M为窗宽，那么，时域相乘相当于频域卷积，频域内，这一窗函数能够分辨出的最近频率肯定不可能小于2*pi/M了，也就是如果数据长度不能满足2*pi/M<|w2-w1|（w2,w1为两个靠的很近的频率），那么在频谱分析时，频谱上将不能分辨出这两个谱，由于w2-w1=2*pi(f2-f1)/fs=2*pi*Δf/fs也就是2*pi/M<2*piΔf/fs，得到Δf的限制为fs/M，这就是窗函数宽度的最小选择，就是说，根据Shannon采样定理确定了采样频率后，要根据靠的最近的谱峰来确定最小的采样长度，这样，所作出来的频谱才能分辨出那两个谱峰，也就是拥有了相应的频率分辨率。
几个例子：
考虑双正弦信号：x = sin(2*pi*10*n)+sin(2*pi*9.8*n);根据Shannon采样定理，采样频率要大于截止频率的两倍，这里选采样频率为80，那么，我们可以看到，Δf为0.2Hz，那么，最小的数据长度为80/0.2=400，但是对正弦信号的频谱分析经验告诉我们，在截断时截断时的数据要包含整周期，并且后面不宜补零以避免频谱泄露（这一点见胡广书《数字信号处理导论》，清华大学出版社），那么，我们要选择至少980个点，才能保含到一个整周期，另外，FFT的经验告诉我们作分析时最好选择2的整数次幂，我们选择靠的最近的1024点。分析结束。
正弦信号根据胡广书老师信号课程的分析，最好大于周期整数4倍或以上.
我觉得在工程上一般取2倍最高频率都不能满足要求了,一般取3-3.5倍比较合适,但是考虑谐频部分,只要取fs>=2*(f1+f2)基本上是不会有问题的....仅个人观点
[CODE]
Fs = 80;
n = 0:1/Fs:1023*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
figure(1);
X = fftshift(fft(x));
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
这是按照我们的分析进行的编程和图形
可以看出这两个谱峰很好的被分辨开来，9.8Hz不在谱线上，所以幅值不为1，以下是一些对比：
[CODE]
Fs = 80;
n = 0:1/Fs:1023*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
figure(1);
subplot(211)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长1024，包含整周期');
n = 0:1/Fs:979*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
subplot(212)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长980，包含整周期');
n = 0:1/Fs:399*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
figure(2);
subplot(211)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长400，不包含整周期');
Fs = 20;
n = 0:1/Fs:1024*1/Fs;
x = sin(2*pi*10*n)+sin(2*pi*9.8*n);
N = length(n);
X = fftshift(fft(x));
subplot(212)
plot((-N/2:N/2-1)*Fs/N,abs(X)*2/N);
grid on;
axis([0 15 0 1]);
title('窗长1024，采样率过小');
喜欢
0
赠金笔
前一篇：fft采样和混叠的概念
后一篇：退耦电容的选择（总结）