1、快速傅立葉變換FFT頻譜分析程序例如，使用ScopFFT快速傅立葉頻譜分析程序對含有噪音的信號（信號成分：振幅6的50Hz正弦波，振幅4的200Hz正弦波，振幅5的250Hz正弦波）進行FFT頻譜分析。 通過ScopFFT快速傅立葉頻譜分析程序進行FFT頻譜分析可以得到下圖的分析結果，可以清晰的分析出振幅約6的50Hz正弦波，振幅約4的200Hz正弦波和振幅約5的250Hz正弦波這三個主要頻率成分。 使用Matlab對采樣數據進行頻譜分析(1)最近做畢設，要對采集到的數據進行頻譜分析。剛開始將所有采樣數據全部送入Matlab進行分析，效果總是很不理想。翻閱課本、搜查網頁，總結出使用Matlab

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

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

4、fft數據采樣個數。一般不指定N，即簡化為Y=fft(b)。Y即為FFT變換后得到的結果，與b的元素數相等，為復數。以頻率為橫坐標，Y數組每個元素的幅值為縱坐標，畫圖即得數據b的幅頻特性；以頻率為橫坐標，Y數組每個元素的角度為縱坐標，畫圖即得數據b的相頻特性。典型頻譜分析M程序舉例如下：clcfs=100;t=0:1/fs:100;N=length(t)-1;%減1使N為偶數%頻率分辨率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);%上面模擬對信號進行采樣，得到采樣數據p，下面對p進行頻譜分析figure(1)plot(t,p);grid ontitle('信號 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('幅值')對于現實中的情況，采樣頻率fs一般都是由采樣儀器決定的，即fs為一個給定的常數；另一方面，為了獲得一定精度的頻譜，對頻率分辨率F有一個人為的規定，一般要求F<0.01，即采樣時間ts>100秒；由采樣時間ts和采樣頻率fs即可決定采樣數據量，即采樣總點數N=fs*ts。這

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

8、2=56.5，不等于整數，即在頻譜圖中找不到準確的刻度，而只能在第56和57個頻率刻度上分開顯示其幅值，這自然就不準確了。因此，請大家在頻譜分析時一定要使F能夠被頻率精度整除。如要求頻率精確度為0.01，則F最大為0.01，也可取值為0.02、0.05、0.001等數據，使0.01/F=整數。而F僅僅由采樣時間ts（也稱數據長度）決定，因此一定要選擇好ts，且要首先確定ts的值。作為驗證，對上面的程序做一個修改：將t=0:1/fs:100;改為t=0:1/fs:83;即ts由100改為83，則F=1/ts由0.01變為0.012。二者分別作出頻譜圖對比如下：上圖1 頻譜圖：ts=100s，F=

9、1/ts=0.01上圖2 頻譜圖：ts=83s，F=1/ts=0.012對比上面兩個圖即可發現，圖2中由于f/F不是整數，在橫坐標中找不到對應的刻度，從而使得各個頻率的幅值泄漏到了其他頻率。總結上面的結論，在保證采樣定理所要求的二倍頻的前提下，并不是采樣頻率fs或采樣點數N越大越好，而是要控制好數據長度ts，使頻率分辨率F滿足頻率精度。FFT實踐及頻譜分析(2)6c995d475575fcd51da4be6.html%      FFT實踐及頻譜分析       &

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

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

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

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

14、s(t,2);x=x(1:99);figure(2);subplot(231);plot(t(1:99),x);%作矩形波的時域波形xlabel('t');ylabel('y');title('矩形波時域波形');grid;%進行FFT變換并做頻譜圖y=fft(x);%進行fft變換mag=abs(y);%求幅值f=(0:length(y)-1)'*fs/length(y);%進行對應的頻率轉換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、%求對數譜ln=log(sq);figure(2);subplot(235);plot(f,ln);xlabel('頻率(Hz)');ylabel('對數譜');title('矩形波對數譜');grid;%用IFFT恢復原始信號xifft=ifft(y);magx=real(xifft);ti=0:length(xifft)-1/fs;figure(2);subplot(236);plot(ti,magx);xlabel('t');ylabel('y');title('通過IFFT轉換的矩形波波形'

17、);grid;%*3.白噪聲*%fs=10;%設定采樣頻率t=-5:0.1:5;x=zeros(1,100);x(50)=100000;figure(3);subplot(231);plot(t(1:100),x);%作白噪聲的時域波形xlabel('t');ylabel('y');title('白噪聲時域波形');grid;%進行FFT變換并做頻譜圖y=fft(x);%進行fft變換mag=abs(y);%求幅值f=(0:length(y)-1)'*fs/length(y);%進行對應的頻率轉換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;%求對數譜ln=log(sq);figure(3);subplot(235);plot(f,ln);xlabel('頻率(Hz)');ylabel('對數譜');title('白噪聲對數譜');grid;%用IFFT恢復原始信號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('通過IFFT轉換的白噪聲波形');grid;FFT的計算結果，對應的頻點計算公式(3)3967a2132a97ddd8da.html1. FFT的頻率分辨率計算公式為：f=fs/N；其中fs為采樣頻率；N為FFT變換的點數2. FFT的計算結果，對應的頻點計算公式：1*fs/N, 2f*s/N, 3f*s/N, N*fs/N。舉例說明：用1KHZ的采樣率采樣信號128點，則FFT結果的128個數據即對應的頻率點分別是1k/128，2k/128, 3k/128, 128k/128HZ。FFT結果的物理意義（轉圈圈）一個FI

21、R濾波器的matlab實現690c33fa19.html問題: 設信號x(t)=sin(2*pi*80*t)+2*sin(2*pi*150*t), 由于某一原因，原始信號被白噪聲污染,實際獲得的信號為x2(t)=x+rand(size(x)，要求設計一個FIR濾波器恢復出原始信號。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、原始信號x');subplot(312);plot(n2,x(1536+n2);title('在原始信號上加上噪聲信號');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('經過濾波器后的信號z');figure(2);subpl

23、ot(311)plot(n1*df,abs(E);title('原始信號頻譜')subplot(312)plot(n1*df,Y);title('噪聲信號的頻譜')subplot(313);plot(n1*df,zk);title('經過濾波器后的信號的頻譜') 濾波后在80150Hz之間噪聲是和濾波器的階數有關，當M越大(程序中用M表示濾波器的階數，濾波器頻響可用freqz來觀看)，80150Hz之間的噪聲越小。用作實驗(兮阿悠1)flag1=input('輸入信號序號:'); 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('輸入信號序號:');  &

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('輸入信號序號:');     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('輸入信號序號:');     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('輸入信號序號:');     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('輸入信號序號:');     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('輸入信號序號:');     N=input('N=&#

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

38、0;       %采樣頻率N = 1024;           %采樣點數t = (0:N-1)/fs;     %采樣時間序列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); %雙邊復數譜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.時域的持續時間范圍應較大；%    2.頻率w與序列k的對應關系：w = 2 * pi* k * df;%    3.采樣頻率 fs 應大于 w/(2*pi) 的2倍%    4.結果曲線的峰值的橫坐標對應的就是 w 和 -w 值fs = 10;       %采樣頻率N = 1024;           %采樣點數t = (0:N-1)/fs;