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

谢谢咯,找了好久了
页: [1] 2 3
查看完整版本: 不同采样率及采样长度的 FFT 分析比较