快速傅立葉變換FFT頻譜分析程序

1、快速傅立葉變換FFT頻譜分析程序例如，使用ScopFFT快速傅立葉頻譜分析程序?qū)性胍舻男盘?hào)（信號(hào)成分：振幅6的50Hz正弦波，振幅4的200Hz正弦波，振幅5的250Hz正弦波）進(jìn)行FFT頻譜分析。 通過(guò)ScopFFT快速傅立葉頻譜分析程序進(jìn)行FFT頻譜分析可以得到下圖的分析結(jié)果，可以清晰的分析出振幅約6的50Hz正弦波，振幅約4的200Hz正弦波和振幅約5的250Hz正弦波這三個(gè)主要頻率成分。 使用Matlab對(duì)采樣數(shù)據(jù)進(jìn)行頻譜分析(1)最近做畢設(shè)，要對(duì)采集到的數(shù)據(jù)進(jìn)行頻譜分析。剛開(kāi)始將所有采樣數(shù)據(jù)全部送入Matlab進(jìn)行分析，效果總是很不理想。翻閱課本、搜查網(wǎng)頁(yè)，總結(jié)出使用Matlab

2、對(duì)采樣數(shù)據(jù)進(jìn)行頻譜分析的理論和方法，特整理于此，和大家分享。1、采樣數(shù)據(jù)導(dǎo)入Matlab采樣數(shù)據(jù)的導(dǎo)入至少有三種方法。第一就是手動(dòng)將數(shù)據(jù)整理成Matlab支持的格式，這種方法僅適用于數(shù)據(jù)量比較小的采樣。第二種方法是使用Matlab的可視化交互操作，具體操作步驟為：File -> Import Data，然后在彈出的對(duì)話框中找到保存采樣數(shù)據(jù)的文件，根據(jù)提示一步一步即可將數(shù)據(jù)導(dǎo)入。這種方法適合于數(shù)據(jù)量較大，但又不是太大的數(shù)據(jù)。據(jù)本人經(jīng)驗(yàn)，當(dāng)數(shù)據(jù)大于15萬(wàn)對(duì)之后，讀入速度就會(huì)顯著變慢，出現(xiàn)假死而失敗。第三種方法，使用文件讀入命令。數(shù)據(jù)文件讀入命令有textread、fscanf、load等，如

3、果采樣數(shù)據(jù)保存在txt文件中，則推薦使用 textread命令。如 a,b=textread('data.txt','%f%*f%f'); 這條命令將data.txt中保存的數(shù)據(jù)三個(gè)三個(gè)分組，將每組的第一個(gè)數(shù)據(jù)送給列向量a，第三個(gè)數(shù)送給列向量b，第二個(gè)數(shù)據(jù)丟棄。命令類(lèi)似于C語(yǔ)言，詳細(xì)可查看其幫助文件。文件讀入命令錄入采樣數(shù)據(jù)可以處理任意大小的數(shù)據(jù)量，且錄入速度相當(dāng)快，一百多萬(wàn)的數(shù)據(jù)不到20秒即可錄入。強(qiáng)烈推薦！2、對(duì)采樣數(shù)據(jù)進(jìn)行頻譜分析頻譜分析自然要使用快速傅里葉變換FFT了，對(duì)應(yīng)的命令即 fft ，簡(jiǎn)單使用方法為：Y=fft(b,N)，其中b即是采樣數(shù)據(jù)，N為

4、fft數(shù)據(jù)采樣個(gè)數(shù)。一般不指定N，即簡(jiǎn)化為Y=fft(b)。Y即為FFT變換后得到的結(jié)果，與b的元素?cái)?shù)相等，為復(fù)數(shù)。以頻率為橫坐標(biāo)，Y數(shù)組每個(gè)元素的幅值為縱坐標(biāo)，畫(huà)圖即得數(shù)據(jù)b的幅頻特性；以頻率為橫坐標(biāo)，Y數(shù)組每個(gè)元素的角度為縱坐標(biāo)，畫(huà)圖即得數(shù)據(jù)b的相頻特性。典型頻譜分析M程序舉例如下：clcfs=100;t=0:1/fs:100;N=length(t)-1;%減1使N為偶數(shù)%頻率分辨率F=1/t=fs/Np=1.3*sin(0.48*2*pi*t)+2.1*sin(0.52*2*pi*t)+1.1*sin(0.53*2*pi*t).+0.5*sin(1.8*2*pi*t)+0.9*sin(2

5、.2*2*pi*t);%上面模擬對(duì)信號(hào)進(jìn)行采樣，得到采樣數(shù)據(jù)p，下面對(duì)p進(jìn)行頻譜分析figure(1)plot(t,p);grid ontitle('信號(hào) p(t)');xlabel('t')ylabel('p')Y=fft(p);magY=abs(Y(1:1:N/2)*2/N;f=(0:N/2-1)'*fs/N;figure(2)%plot(f,magY);h=stem(f,magY,'fill','-');set(h,'MarkerEdgeColor','red',

6、9;Marker','*')grid ontitle('頻譜圖 （理想值：0.48Hz,1.3、0.52Hz,2.1、0.53Hz,1.1、1.8Hz,0.5、2.2Hz,0.9） ');xlabel('f (Hz)')ylabel('幅值')對(duì)于現(xiàn)實(shí)中的情況，采樣頻率fs一般都是由采樣儀器決定的，即fs為一個(gè)給定的常數(shù)；另一方面，為了獲得一定精度的頻譜，對(duì)頻率分辨率F有一個(gè)人為的規(guī)定，一般要求F<0.01，即采樣時(shí)間ts>100秒；由采樣時(shí)間ts和采樣頻率fs即可決定采樣數(shù)據(jù)量，即采樣總點(diǎn)數(shù)N=fs*ts。這

7、就從理論上對(duì)采樣時(shí)間ts和采樣總點(diǎn)數(shù)N提出了要求，以保證頻譜分析的精準(zhǔn)度。3、數(shù)據(jù)長(zhǎng)度的選擇頻率分辨率F，顧名思義就是頻譜中能夠區(qū)分出的最小頻率刻度。如F=0.01，則頻譜圖中橫坐標(biāo)頻率的最小刻度為0.01，即0.02Hz和 0.03Hz是沒(méi)有準(zhǔn)確數(shù)據(jù)的，但Matlab在畫(huà)圖時(shí)對(duì)其進(jìn)行了插值，故而plot作圖時(shí)看到的頻譜是連續(xù)的。但用stem來(lái)作圖就可以看出頻率是離散的，stem對(duì)了解F的含義非常有幫助。由此，我們可以進(jìn)一步思考。如果信號(hào)所包含的頻率分量不是F的整數(shù)倍，那么這個(gè)頻率分量就不會(huì)得到正確的反映。如信號(hào)包含1.13Hz頻率分量，而 F=1/ts=fs/N=0.02，則1.13/0.0

8、2=56.5，不等于整數(shù)，即在頻譜圖中找不到準(zhǔn)確的刻度，而只能在第56和57個(gè)頻率刻度上分開(kāi)顯示其幅值，這自然就不準(zhǔn)確了。因此，請(qǐng)大家在頻譜分析時(shí)一定要使F能夠被頻率精度整除。如要求頻率精確度為0.01，則F最大為0.01，也可取值為0.02、0.05、0.001等數(shù)據(jù)，使0.01/F=整數(shù)。而F僅僅由采樣時(shí)間ts（也稱(chēng)數(shù)據(jù)長(zhǎng)度）決定，因此一定要選擇好ts，且要首先確定ts的值。作為驗(yàn)證，對(duì)上面的程序做一個(gè)修改：將t=0:1/fs:100;改為t=0:1/fs:83;即ts由100改為83，則F=1/ts由0.01變?yōu)?.012。二者分別作出頻譜圖對(duì)比如下：上圖1 頻譜圖：ts=100s，F(xiàn)=

9、1/ts=0.01上圖2 頻譜圖：ts=83s，F(xiàn)=1/ts=0.012對(duì)比上面兩個(gè)圖即可發(fā)現(xiàn)，圖2中由于f/F不是整數(shù)，在橫坐標(biāo)中找不到對(duì)應(yīng)的刻度，從而使得各個(gè)頻率的幅值泄漏到了其他頻率?？偨Y(jié)上面的結(jié)論，在保證采樣定理所要求的二倍頻的前提下，并不是采樣頻率fs或采樣點(diǎn)數(shù)N越大越好，而是要控制好數(shù)據(jù)長(zhǎng)度ts，使頻率分辨率F滿足頻率精度。FFT實(shí)踐及頻譜分析(2)6c995d475575fcd51da4be6.html%      FFT實(shí)踐及頻譜分析       &

10、#160;                  %*%*1.正弦波*%fs=100;%設(shè)定采樣頻率N=128;n=0:N-1;t=n/fs;f0=10;%設(shè)定正弦信號(hào)頻率%生成正弦信號(hào)x=sin(2*pi*f0*t);figure(1);subplot(231);plot(t,x);%作正弦信號(hào)的時(shí)域波形xlabel('t');ylabel('y');title('正弦信號(hào)y=2*pi*1

11、0t時(shí)域波形');grid;%進(jìn)行FFT變換并做頻譜圖y=fft(x,N);%進(jìn)行fft變換mag=abs(y);%求幅值f=(0:length(y)-1)'*fs/length(y);%進(jìn)行對(duì)應(yīng)的頻率轉(zhuǎn)換figure(1);subplot(232);plot(f,mag);%做頻譜圖axis(0,100,0,80);xlabel('頻率(Hz)');ylabel('幅值');title('正弦信號(hào)y=2*pi*10t幅頻譜圖N=128');grid;%求均方根譜sq=abs(y);figure(1);subplot(233);p

12、lot(f,sq);xlabel('頻率(Hz)');ylabel('均方根譜');title('正弦信號(hào)y=2*pi*10t均方根譜');grid;%求功率譜power=sq.2;figure(1);subplot(234);plot(f,power);xlabel('頻率(Hz)');ylabel('功率譜');title('正弦信號(hào)y=2*pi*10t功率譜');grid;%求對(duì)數(shù)譜ln=log(sq);figure(1);subplot(235);plot(f,ln);xlabel('

13、;頻率(Hz)');ylabel('對(duì)數(shù)譜');title('正弦信號(hào)y=2*pi*10t對(duì)數(shù)譜');grid;%用IFFT恢復(fù)原始信號(hào)xifft=ifft(y);magx=real(xifft);ti=0:length(xifft)-1/fs;figure(1);subplot(236);plot(ti,magx);xlabel('t');ylabel('y');title('通過(guò)IFFT轉(zhuǎn)換的正弦信號(hào)波形');grid;%*2.矩形波*%fs=10;%設(shè)定采樣頻率t=-5:0.1:5;x=rectpul

14、s(t,2);x=x(1:99);figure(2);subplot(231);plot(t(1:99),x);%作矩形波的時(shí)域波形xlabel('t');ylabel('y');title('矩形波時(shí)域波形');grid;%進(jìn)行FFT變換并做頻譜圖y=fft(x);%進(jìn)行fft變換mag=abs(y);%求幅值f=(0:length(y)-1)'*fs/length(y);%進(jìn)行對(duì)應(yīng)的頻率轉(zhuǎn)換figure(2);subplot(232);plot(f,mag);%做頻譜圖xlabel('頻率(Hz)');ylabel(&

15、#39;幅值');title('矩形波幅頻譜圖');grid;%求均方根譜sq=abs(y);figure(2);subplot(233);plot(f,sq);xlabel('頻率(Hz)');ylabel('均方根譜');title('矩形波均方根譜');grid;%求功率譜power=sq.2;figure(2);subplot(234);plot(f,power);xlabel('頻率(Hz)');ylabel('功率譜');title('矩形波功率譜');grid;

16、%求對(duì)數(shù)譜ln=log(sq);figure(2);subplot(235);plot(f,ln);xlabel('頻率(Hz)');ylabel('對(duì)數(shù)譜');title('矩形波對(duì)數(shù)譜');grid;%用IFFT恢復(fù)原始信號(hào)xifft=ifft(y);magx=real(xifft);ti=0:length(xifft)-1/fs;figure(2);subplot(236);plot(ti,magx);xlabel('t');ylabel('y');title('通過(guò)IFFT轉(zhuǎn)換的矩形波波形'

17、);grid;%*3.白噪聲*%fs=10;%設(shè)定采樣頻率t=-5:0.1:5;x=zeros(1,100);x(50)=100000;figure(3);subplot(231);plot(t(1:100),x);%作白噪聲的時(shí)域波形xlabel('t');ylabel('y');title('白噪聲時(shí)域波形');grid;%進(jìn)行FFT變換并做頻譜圖y=fft(x);%進(jìn)行fft變換mag=abs(y);%求幅值f=(0:length(y)-1)'*fs/length(y);%進(jìn)行對(duì)應(yīng)的頻率轉(zhuǎn)換figure(3);subplot(232

18、);plot(f,mag);%做頻譜圖xlabel('頻率(Hz)');ylabel('幅值');title('白噪聲幅頻譜圖');grid;%求均方根譜sq=abs(y);figure(3);subplot(233);plot(f,sq);xlabel('頻率(Hz)');ylabel('均方根譜');title('白噪聲均方根譜');grid;%求功率譜power=sq.2;figure(3);subplot(234);plot(f,power);xlabel('頻率(Hz)')

19、;ylabel('功率譜');title('白噪聲功率譜');grid;%求對(duì)數(shù)譜ln=log(sq);figure(3);subplot(235);plot(f,ln);xlabel('頻率(Hz)');ylabel('對(duì)數(shù)譜');title('白噪聲對(duì)數(shù)譜');grid;%用IFFT恢復(fù)原始信號(hào)xifft=ifft(y);magx=real(xifft);ti=0:length(xifft)-1/fs;figure(3);subplot(236);plot(ti,magx);xlabel('t'

20、);ylabel('y');title('通過(guò)IFFT轉(zhuǎn)換的白噪聲波形');grid;FFT的計(jì)算結(jié)果，對(duì)應(yīng)的頻點(diǎn)計(jì)算公式(3)3967a2132a97ddd8da.html1. FFT的頻率分辨率計(jì)算公式為：f=fs/N；其中fs為采樣頻率；N為FFT變換的點(diǎn)數(shù)2. FFT的計(jì)算結(jié)果，對(duì)應(yīng)的頻點(diǎn)計(jì)算公式：1*fs/N, 2f*s/N, 3f*s/N, N*fs/N。舉例說(shuō)明：用1KHZ的采樣率采樣信號(hào)128點(diǎn)，則FFT結(jié)果的128個(gè)數(shù)據(jù)即對(duì)應(yīng)的頻率點(diǎn)分別是1k/128，2k/128, 3k/128, 128k/128HZ。FFT結(jié)果的物理意義（轉(zhuǎn)圈圈）一個(gè)FI

21、R濾波器的matlab實(shí)現(xiàn)690c33fa19.html問(wèn)題: 設(shè)信號(hào)x(t)=sin(2*pi*80*t)+2*sin(2*pi*150*t), 由于某一原因，原始信號(hào)被白噪聲污染,實(shí)際獲得的信號(hào)為x2(t)=x+rand(size(x)，要求設(shè)計(jì)一個(gè)FIR濾波器恢復(fù)出原始信號(hào)。t=0.0:0.001:2.047;x=sin(2*pi*80.*t)+2*sin(2*pi*150.*t);x1=x+randn(size(x);l=length(x);N=1:l;n1=1:1024;n2=1:200;M=64;subplot(311);plot(n2,x(1536+n2);title('

22、原始信號(hào)x');subplot(312);plot(n2,x(1536+n2);title('在原始信號(hào)上加上噪聲信號(hào)');Y=fft(x1);E=fft(x);E=abs(E(n1);Y=abs(Y(n1);fs=1000;df=fs/l;Wn=75 85 145 155/500;b=fir1(M,Wn);%freqz(b,1,512);z=filter(b,1,x1); %zk=fft(z);zk=abs(zk(n1);subplot(313)plot(n2,z(1536+n2);title('經(jīng)過(guò)濾波器后的信號(hào)z');figure(2);subpl

23、ot(311)plot(n1*df,abs(E);title('原始信號(hào)頻譜')subplot(312)plot(n1*df,Y);title('噪聲信號(hào)的頻譜')subplot(313);plot(n1*df,zk);title('經(jīng)過(guò)濾波器后的信號(hào)的頻譜') 濾波后在80150Hz之間噪聲是和濾波器的階數(shù)有關(guān)，當(dāng)M越大(程序中用M表示濾波器的階數(shù)，濾波器頻響可用freqz來(lái)觀看)，80150Hz之間的噪聲越小。用作實(shí)驗(yàn)(兮阿悠1)flag1=input('輸入信號(hào)序號(hào):'); N=input('N='); wh

24、ile flag1=1     n=1:1:N     x=1,1,1,1;     X=(x,N);     figure(1);     subplot(2,1,1);     stem(x,'r');     title('x(n)的圖形');     ylabel('x(n)');  

25、0;  xlabel('n');     subplot(2,1,2);     stem(abs(X),'r');     title('|X(k)|的圖形');     ylabel('|X(k)|');     xlabel('k');     flag1=input('輸入信號(hào)序號(hào):');  &

26、#160;  N=input('N='); end while flag1=2     n=1:8     x=1,2,3,4,4,3,2,1     X=(x,N);     figure(2);     subplot(2,1,1);     stem(x,'r');     title('x(n)的圖形');  &

27、#160;  ylabel('x(n)');     xlabel('n');     subplot(2,1,2);     stem(abs(X),'r');     title('|X(k)|的圖形');     ylabel('|X(k)|');     xlabel('k');   

28、  flag1=input('輸入信號(hào)序號(hào):');     N=input('N='); end while flag1=3     n=1:8     x=4,3,2,1,1,2,3,4     X=(x,N);     figure(3);     subplot(2,1,1);     stem(x,'r');  

29、   title('x(n)的圖形');     ylabel('x(n)');     xlabel('n');     subplot(2,1,2);     stem(abs(X),'r');     title('|X(k)|的圖形');     ylabel('|X(k)|');  

30、   xlabel('k');     flag1=input('輸入信號(hào)序號(hào):');     N=input('N='); end while flag1=4     n=1:N     x=cos(n*pi/4)     X=(x,N);     figure(4);     subplot(2,1,1); 

31、0;   stem(x,'r');     title('x(n)的圖形');     ylabel('x(n)');     xlabel('n');     subplot(2,1,2);     stem(abs(X),'r');     title('|X(k)|的圖形');  

32、60;  ylabel('|X(k)|');     xlabel('k');     flag1=input('輸入信號(hào)序號(hào):');     N=input('N='); end while flag1=5      n=1:N      x=sin(n*pi/8);      X=(x,N);  &#

33、160;   figure(5);     subplot(2,1,1);     stem(x,'r');     title('x(n)的圖形');     ylabel('x(n)');     xlabel('n');     subplot(2,1,2);     stem(abs(X),

34、9;r');     title('|X(k)|的圖形');     ylabel('|X(k)|');     xlabel('k');     flag1=input('輸入信號(hào)序號(hào):');     N=input('N='); end while flag1=6     n=1:N     x=

35、cos(pi*n/8)+cos(pi*n/4)+cos(5*pi*n/16);     X=(x,N);     figure(6);     subplot(2,1,1);     stem(x,'r');     title('x(n)的圖形');     ylabel('x(n)');     xlabel('n');

36、     subplot(2,1,2);     stem(abs(X),'r');     title('|X(k)|的圖形');     ylabel('|X(k)|');     xlabel('k');     flag1=input('輸入信號(hào)序號(hào):');     N=input('N=&#

37、39;); end關(guān)于FFT的頻譜對(duì)應(yīng)關(guān)系(兮阿悠2)0f731add187.html根據(jù)論壇上面的帖子重新總結(jié)了一下，這下應(yīng)該完整了。有兩種方案，均可成功顯示調(diào)用FFt后的頻譜圖（主要是突出頻譜圖橫坐標(biāo)和原信號(hào)的一致性）% 方案1：“x = a*cos(2*pi*w*t)”的形式：% -% 注意：1.時(shí)域的持續(xù)時(shí)間范圍應(yīng)較大；%    2.頻率w與序列k的對(duì)應(yīng)關(guān)系：w =   k * df;%    3.采樣頻率 fs 應(yīng)大于 w 的2倍%    4.結(jié)果曲線的峰值的橫坐標(biāo)對(duì)應(yīng)的就是 w 和 -w 值fs = 1

38、0;       %采樣頻率N = 1024;           %采樣點(diǎn)數(shù)t = (0:N-1)/fs;     %采樣時(shí)間序列sa = 0.75;w = 4;x = a*cos(2*pi*w*t);subplot(2,1,1);plot(t, x);xlabel('t/s');xf = fft(x,N)/N; xf = fftshift(xf); %雙邊復(fù)數(shù)譜df = fs/N;       %

39、頻率分辨率Hzf = (-N/2+1:N/2)*df; %頻域序列subplot(2,1,2);plot(f, abs(xf);xlabel('f/Hz');% 方案2：“x = a*cos(w*t)”的形式：-% 注意：1.時(shí)域的持續(xù)時(shí)間范圍應(yīng)較大；%    2.頻率w與序列k的對(duì)應(yīng)關(guān)系：w = 2 * pi* k * df;%    3.采樣頻率 fs 應(yīng)大于 w/(2*pi) 的2倍%    4.結(jié)果曲線的峰值的橫坐標(biāo)對(duì)應(yīng)的就是 w 和 -w 值fs = 10;       %采樣頻率N = 1024;           %采樣點(diǎn)數(shù)t = (0:N-1)/fs;