论坛 › 信号处理方法 › 不同采样率及采样长度的 FFT 分析比较

xuefei01 发表于 2006-12-28 21:38

不同采样率及采样长度的 FFT 分析比较

不同采样率及采样长度的 FFT 分析比较

[ 本帖最后由 zhlong 于 2007-6-4 17:51 编辑 ]

shenyongjun 发表于 2006-12-29 10:27

嗯，资料不错。如果现在正在学习《现代谱估计》之类的课程的话，这可以作为一个作业进行练习！

xuefei01 发表于 2007-1-1 12:35

补程序

clear
clf
f1=1000;% f1=1kHz
f2=2500;% f2=2.5kHz
f3=3000;% f3=3kHz
fs=10000; % 采样频率f0=10kHZ
T=1/fs;   % 采样周期
t=(0:100).*T;
y=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);
plot(t,y);
N=20;   % 采样点数
Nfft=20;% FFT 点数不足后面补零
W0=fft(y(1:N),Nfft);
Aw=2*abs(W0)/N;
df=fs/Nfft;% [0,fs） 频点间隔

%改变采样点数 和fft变换点数 比较分析
N1=40;   % 采样点数
Nfft1=128;% FFT 点数不足后面补零
W1=fft(y(1:N1),Nfft1);
Aw1=2*abs(W1)/N1;
df1=fs/Nfft1;% [0,fs） 频点间隔

figure(2);
plot((0:Nfft-1)*df,Aw,'-r');
hold on
plot((0:Nfft1-1)*df1,Aw1);

csujds 发表于 2007-1-3 19:25

这行是什么意思，为什么×2/N

Aw=2*abs(W0)/N;

heliping1981 发表于 2007-1-5 11:10

这么好的程序,谁解释一下撒!

kiddjiang 发表于 2007-1-6 18:52

乘以2是变换成单边谱
除以N是得到各频率成分的正确系数
否则的话，幅度会随着N变化而变化

shirley329 发表于 2007-1-7 06:29

x(n)=cos(0.48*pi*n)+cos(0.52*pi*n),利用MATLAB程序求如下X（exp(jw)）,X(k)。
(1) 取x(n)的前10点数据，求N=10点的X（exp(jw)）,X(k)并作图。
(2) 取（1）中的x(n)补零至100点，求N=100点的X（exp(jw)）,X(k)并作图。
(3) 取x(n)的前100点数据，求N=100点的X（exp(jw)）,X(k)并作图。
(4) 取x(n)的前128点数据，求N=128点的X（exp(jw)）,X(k)并作图。
(5) 取x(n)的前50点数据，求N=100点的X（exp(jw)）,X(k)并作图。
(6) 讨论以上五种情况的区别。

2．(1) matlab程序如下：
figure(1) n=;
x=cos(0.48*pi*n)+cos(0.52*pi*n);
w=*2*pi/500;%0-2*pi区域分为501点
subplot(3,1,1);
stem(n,x);
title('x(n) (0<=n<=9');
xlabel('n');
axis();% axis()
line(,);
x1=fft(x);
magx1=abs(x1(1:1:10));
x=x*exp(-j*n'*w);%内部的矩阵维数必须一致
magx=abs(x);
k1=0:1:9;
w1=2*pi/10*k1;
subplot(3,1,2);
plot(w/pi,magx);
title('DTFT');
xlabel('w');
axis();
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis();

图（1） （2）


程序如下：
n=;
y=cos(0.48*pi*n)+cos(0.52*pi*n);
n1=;
x=;
subplot(3,1,1);
stem(n1,x);
title('x(n) (0<=n<=9+90zeros');
xlabel('n');
axis();% axis()
line(,);
w=*2*pi/500;%0-2*pi区域分为501点
x1=fft(x);
magx1=abs(x1(1:1:51));
x=x*exp(-j*n1'*w);
magx=abs(x);
k1=;
w1=2*pi/100*k1;
subplot(3,1,2);
plot(w/pi,magx);
title('DTFT');
xlabel('w');
axis();
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis();
结果如下： 图（2） （3）


取x(n)的前100点数据，程序如下：
n=;
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(3,1,1);
stem(n,x);
title('x(n) (0<=n<=99');
xlabel('n');
axis();%axis()
line(,);
w=*2*pi/500;%0-2*pi区域分为501点
x1=fft(x);
magx1=abs(x1(1:1:50));
x=x*exp(-j*n'*w);
magx=abs(x);
k1=0:1:49;
w1=2*pi/100*k1;
subplot(3,1,2);
plot(w/pi,magx);
title('DTFT');
xlabel('w');
axis();
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis();
结果如下： 图（3） （4）


n=;
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(3,1,1); stem(n,x);
title('x(n) (0<=n<=127');
xlabel('n');
axis();% axis()
line(,);
w=*2*pi/500;%0-2*pi区域分为501点
x1=fft(x);
magx1=abs(x1(1:1:64));
x=x*exp(-j*n'*w);
magx=abs(x); k1=0:1:63;
w1=2*pi/128*k1;
subplot(3,1,2);
plot(w/pi,magx);
title('DTFT');
xlabel('w');
axis();
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis();
图（4） （5）


n=; x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(3,1,1);
stem(n,x);
title('x(n) (0<=n<=49');
xlabel('n');
axis();% axis()
line(,);
w=*2*pi/500;%0-2*pi区域分为501点
x1=fft(x);
magx1=abs(x1(1:1:50));
x=x*exp(-j*n'*w);
magx=abs(x);
k1=0:1:49;
w1=2*pi/50*k1;
subplot(3,1,2);
plot(w/pi,magx);
title('DTFT');
xlabel('w');
axis();
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis();
图（5）


2． 分析：由图（1）可见，由于截断函数的频谱混叠作用，X（K）不能正确分辨w1=0.48*pi,w2=0.52*pi这两个频率分量。
由图（2）可见，虽然x(n)补零至100点，X（K）的密度，截断函数的频谱混叠作用没有改变，这时的物理分辨率使X(K)仍不能正确分辨w1=0.48*pi,w2=0.52*pi这两个频率分量。
由图（3）可见，截断函数的加宽且为周期序列的整数倍，改变了频谱混叠作用，提高了“物理”分辨率使X（K）能正确分辨w1=0.48*pi,w2=0.52*pi这两个频率分量。
由图（4）可见，截断函数虽进一步加宽，但不是周期序列的整数倍，所以尽管 X（K）能正确分辨w1=0.48*pi,w2=0.52*pi这两个频率分量，但还是呈现除了频谱泄露。
由图（5）可见，截断函数的宽度正好为序列的周期，即这时的“物理”分辨率使X（K）正好能正确分辨w1=0.48*pi,w2=0.52*pi这两个频率分量。

物理分辨率低与频谱的混叠有关，而频谱的混叠正是由截断造成的。若由两个不同频率的周期余弦信号组合的x(n),仍是一个周期余弦信号的话，其新周期N0=2*pi/w’,w’=|w2-w1|,当N0>Np时，必有信息损失，导致频谱混叠，严重的就无法分辨原有谱峰，本例中周期N=50,所以n’>=50. 因而图（1），（2）中都出现了频谱混叠。要改变频谱分辨率，就必须加宽截断函数的时宽Np,图（3），（4），（5）都做到了这一点，因而都有效的克服了频谱混叠效应。在数据长度Np一定的情况下，尽管可在x(n)后面加零，改变频谱的频率取样间隔，但改变的仅是频谱的频率取样密度，而无法改变频率分辨率，这一点从图（2）中可以看出来。
由此，我们引出另外一个问题：栅栏效应。序列x(n)的DTFT频谱是连续的（当然，同时它也是周期的），而DFT所求的是这个连续谱上的均匀抽样。如果用X(k)去近似，就一定意义上来讲，就好象是在栅栏的一边通过栅栏的缝隙(对应离散点)去观看另一边的景象(对应连续频谱)，只能在离散点的地方看到真实的景象，因此，那些被栅栏挡住的（频谱）部分是看不到的，这就有可能漏掉一些较大频率分量。我们称这种现象为"栅栏效应"。 当然，在实际问题中，"大的频谱分量"被挡住的情形还是很少的，栅栏效应并不是一个很严重的问题。尽管如此，我们还是有必要讨论清楚如何避免或者说减少这种栅栏效应。　　
那么，到底怎么才能减少这种栅栏效应呢？ 　　
知道了栅栏产生原因，我们才能采取有针对性的方法来减少它。从根本上讲，用离散的DFT谱来近似连续的DTFT谱，误差总是有的，即从理论上，栅栏效应是不可能消除的。不过，也正是因为DFT离散谱是DTFT连续谱的抽样，为了提高近似程度-，也即减少栅栏效应，我们可以在DTFT连续谱上多抽取些频率样本。这样，谱线变密了，那些原先看不到的频谱分量，就有可能看到了。　　
那么，怎样才能实现在DTFT连续谱上抽取更多的频率点呢？根据DFT与DTFT的关系，只要增加下列DFT的计算式中的N值即可 考虑到增加N值，会使得L=N0，但Np!=MN0时，可以看到频谱泄漏,如图（4）。
所谓频谱泄漏,就是信号频谱中各谱线之间相互影响!使测量结果偏离实际值!同时在谱线两侧其他频率点上出现一些幅值较小的假谱.简单说来!造成频谱泄漏的原因是采样频率与信号频率不同步!造成周期采样信号的相位在始端和终端不连续. 对于频谱泄漏，只要保证窗口函数的宽度为基波周期的整数倍，就可以避免泄漏效应的产生。
其解决办法有二，一是采用适当的窗函数来降低泄漏效应的影响，但是，这种方法同时也增加了计算量。对于大数据量的数据处理而言是不合适的；其二，也是最实用、最有效的解决办法，设计有效的频率跟踪电路，使采样频率实时跟踪信号的基波频率。也就是根据采样时的基波频率来确定采样间隔，从而从根本上解决频谱泄漏效应。

http://i.ttsite.com/u/Mon_0701/Day_2/467007_18975_52e227c04a64c2b.jpg
http://i.ttsite.com/u/Mon_0701/Day_2/467007_18975_39f69afd2b8d7e7.jpg
http://i.ttsite.com/u/Mon_0701/Day_2/467007_18975_cce6cf3a51ae7e3.jpg
http://i.ttsite.com/u/Mon_0701/Day_2/467007_18975_d33dc40ab50eaaf.jpg
http://i.ttsite.com/u/Mon_0701/Day_2/467007_18975_810eec1c321ba46.jpg

http://s329.ttsite.com/read.php?tid=1764099&fid=467007&page=1&toread=1

[ 本帖最后由 shirley329 于 2007-1-7 06:32 编辑 ]

xuefei01 发表于 2007-1-7 13:54

shirley329 分析的相当好

共同学习了

lingchen126 发表于 2007-1-7 21:57

xuefei01是不是广州精信的？

请联系

xuefei01 发表于 2007-1-8 12:05

建议 给shirley329 增加威望值

shirley329 发表的内容不错建议增加威望

hequn 发表于 2007-3-5 14:50

看不到图！

请shirley329重新上传图1～5，谢谢！

hequn 发表于 2007-3-6 15:46

图3没有

请shirley329重新上传图3，谢谢！

piko11 发表于 2010-4-6 19:44

好东西，大家学习。。。

江湖小生 发表于 2010-4-7 16:41

x=x*exp(-j*n'*w);
这里面的运算之前的x是（1，100）的，n'是（100，1）的，w是（501，1）的，满足内部的矩阵维数必须一致 的话，就应该前面两个相乘以后再和后面的向量运算，但是习惯是括号里面的先相乘的啊，这样就不满足了。本人是新手，求高手指导

sail_fan 发表于 2010-5-12 14:20

谢谢咯，找了好久了

查看完整版本: 不同采样率及采样长度的 FFT 分析比较