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

xuefei01

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

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

shenyongjun

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

xuefei01

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

Aw=2*abs(W0)/N;

heliping1981

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

kiddjiang

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

shirley329

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=[0:1:9];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
w=[0:1:500]*2*pi/500;%0-2*pi区域分为501点
subplot(3,1,1);
stem(n,x);
title('x(n) (0<=n<=9');
xlabel('n');
axis([0,10,-2.5,2.5]);% axis([xmin xmax ymin ymax])
line([0,10],[0,0]);
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([0,1,0,10]);
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis([0,1,0,10]);

图（1） （2）


程序如下：
n=[0:1:9];
y=cos(0.48*pi*n)+cos(0.52*pi*n);
n1=[0:1:99];
x=[y(1:1:10),zeros(1,90)];
subplot(3,1,1);
stem(n1,x);
title('x(n) (0<=n<=9+90zeros');
xlabel('n');
axis([0,100,-2.5,2.5]);% axis([xmin xmax ymin ymax])
line([0,100],[0,0]);
w=[0:1:500]*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=[0:1:50];
w1=2*pi/100*k1;
subplot(3,1,2);
plot(w/pi,magx);
title('DTFT');
xlabel('w');
axis([0,1,0,10]);
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis([0,1,0,10]);
结果如下： 图（2） （3）


取x(n)的前100点数据，程序如下：
n=[0:1:99];
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([0,100,-2.5,2.5]);%axis([xmin xmax ymin ymax])
line([0,100],[0,0]);
w=[0:1:500]*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([0,1,0,55]);
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis([0,1,0,55]);
结果如下： 图（3） （4）


n=[0:1:127];
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([0,127,-2.5,2.5]);% axis([xmin xmax ymin ymax])
line([0,100],[0,0]);
w=[0:1:500]*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([0,1,0,70]);
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis([0,1,0,70]);
图（4） （5）


n=[0:1:49]; 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([0,49,-2.5,2.5]);% axis([xmin xmax ymin ymax])
line([0,50],[0,0]);
w=[0:1:500]*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([0,1,0,30]);
subplot(3,1,3);
stem(w1/pi,magx1);
title('X(K)');
xlabel('frequency in pi units');
axis([0,1,0,30]);
图（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://s329.ttsite.com/read.php? ... page=1&toread=1

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

xuefei01

共同学习了

lingchen126

请联系

xuefei01

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

hequn

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

hequn

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

piko11

好东西，大家学习。。。

江湖小生

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

sail_fan

谢谢咯，找了好久了