the sample frenqy and window with and fft

16 篇文章 1 订阅
MATLAB中使用FFT做频谱分析时频率分辨率问题
最近做FFT时,使用的采样频率和信号长度的取舍一直没有搞清楚,后来在论坛上发了一个贴子《总结一下使用FFT和维纳-辛钦定理求解PSD问题》(讨论见http://www.ilovematlab.cn/thread-27150-1-1.html特别感谢会员songzy41,他的问题给了我很大启示),跟帖中给了我不少启示,并且让我对“频率分辨率”这个概念有了更深入的理解。再次一并感谢论坛的高手们。
频率分辨率,顾名思义,就是将信号中两个靠的很近的频谱分开的能力。
信号x(t)长度为Ts,通过傅氏变换后得到X,其频率分辨率为Δf=1/THz),若经过采样后,假设采样频率为fs=1/Ts,而进行频谱分析时要将这个无穷长的序列使用窗函数截断处理,假设使用矩形窗,我们知道,矩形窗的频谱为sinc函数,主瓣宽度可以定义为2*pi/MM为窗宽,那么,时域相乘相当于频域卷积,频域内,这一窗函数能够分辨出的最近频率肯定不可能小于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,那么,我们可以看到,Δf0.2Hz,那么,最小的数据长度为0.2/80=400,但是对正弦信号的频谱分析经验告诉我们,在截断时截断时的数据要包含整周期,并且后面不宜补零以避免频谱泄露(这一点见胡广书《数字信号处理导论》,清华大学出版社),那么,我们要选择至少980个点,才能保含到一个整周期,另外,FFT的经验告诉我们作分析时最好选择2的整数次幂,我们选择靠的最近的1024点。分析结束。

[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,采样率过小');
总结一下FFT和维纳辛钦定理求解PSD的问题-功率谱图.rar 早上在论坛上问了两个问题, 一个是关于FFT求频谱时纵坐标的值问题 https://www.ilovematlab.cn/thread-27092-1-1.html 还有一个是用维纳辛钦定理求解PSD时出现的问题 https://www.ilovematlab.cn/thread-27133-1-1.html 经过达人们的指点,和自己的总结,获得一点心得,在这里与大家分享一下:) 1.FFT求频谱 [CODE] Fs = 40; n = 0:1/Fs:159*1/Fs; x = sin sin; N = length; X = fftshift); Px1 = X.*conj/N; plot*Fs/N,Px1); grid on; axis title; 首先,fftshift的问题,以前上数字信号处理时,老师专门给提出了这个函数,但是我发现论坛里好多不太明白这个函数意义的,OO~,一般,fft得到的是频谱范围在【0-2*pi】范围内的频谱,以高频pi为中心,但是一般使用过程中,使用的频谱习惯以低频0为中心,fftshift的功能就是将频谱进行移位,使之在【-pi,pi】之间; 另外,纵坐标的问题,版主edifier2008提示说用/N的方法归一化,我试了一下,每次采样长度变大时,纵坐标的整体值都会变大,/N之后,值变为1之内了,但是并不是理论算法中得到的1. 图形如下: fft.jpg fft 2.维纳辛钦定理求解功率谱的问题 [CODE] Fs = 40; n = 0:1/Fs:159*1/Fs; x = sin sin; N = length; Rx = xcorr; Px2 = fftshift); plot*Fs/,abs); grid on; axis title; 图形如下: fftwei.jpg 程序中可以看出,也要使用fftshift对fft得到的频谱进行移位以得到以低频0为中心的频谱,另外,得到的功率谱纵轴值特别大,是不是也需要除以采样长度,我试了一下,仍然是很大,个人认为,在MATLAB中计算自相关函数以及计算FFT时,都没有对加和进行归一,将/N这一个系数可能都给省略掉了。 此外,我在很多教材里面看了不少里面的例题,都没有注意纵轴值的问题,我觉得在进行频谱分析,重点在于频率点,以及相近频率点的谱图是不是能够分辨出来,而对于各谱的大小,有个相对的比较即可。 不当之处,还望大家给与指正,:) :victory:
MATLAB中使用FFT做频谱分析时频率分辨率问题-频率分辨率.rar MATLAB中使用FFT做频谱分析时频率分辨率问题 最近做FFT时,使用的采样频率和信号长度的取舍一直没有搞清楚,后来在论坛上发了一个贴子《总结一下使用FFT和维纳-辛钦定理求解PSD问题》(讨论见https://www.ilovematlab.cn/thread-27150-1-1.html,特别感谢会员songzy41,他的问题给了我很大启示),跟帖中给了我不少启示,并且让我对“频率分辨率”这个概念有了更深入的理解。再次一并感谢论坛的高手们。 频率分辨率,顾名思义,就是将信号中两个靠的很近的频谱分开的能力。 信号x长度为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/fs=2*pi*Δf/fs也就是2*pi/M<2*piΔf/fs,得到Δf的限制为fs/M,这就是窗函数宽度的最小选择,就是说,根据Shannon采样定理确定了采样频率后,要根据靠的最近的谱峰来确定最小的采样长度,这样,所作出来的频谱才能分辨出那两个谱峰,也就是拥有了相应的频率分辨率。 几个例子: 考虑双正弦信号:x = sin sin;根据Shannon采样定理,采样频率要大于截止频率的两倍,这里选采样频率为80,那么,我们可以看到,Δf为0.2Hz,那么,最小的数据长度为0.2/80=400,但是对正弦信号的频谱分析经验告诉我们,在截断时截断时的数据要包含整周期,并且后面不宜补零以避免频谱泄露(这一点见胡广书《数字信号处理导论》,清华大学出版社),那么,我们要选择至少980个点,才能保含到一个整周期,另外,FFT的经验告诉我们作分析时最好选择2的整数次幂,我们选择靠的最近的1024点。分析结束。 [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; figure; X = fftshift); plot*Fs/N,abs*2/N); grid on; axis; 这是按照我们的分析进行的编程和图形 zheng.jpg 可以看出这两个谱峰很好的被分辨开来,9.8Hz不在谱线上,所以幅值不为1,以下是一些对比: [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:979*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:399*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; Fs = 20; n = 0:1/Fs:1024*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; 结果如下: 1024.jpg 400.jpg 这是我在做FFT以及论坛中的问题时所得到的一点启发,不当之处还请大家指正。OO~ 频率分辨率.rar 为了方便大家,我将doc版报告和m文件一起上传,和帖子内容一样。OO~
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值