用MATLAB對信號做頻譜分析


1.首先學習下傅里葉變換的東西。學高數的時候老師只是將傅里葉變換簡單的說了下,並沒有深入的講解。而現在看來,傅里葉變換似乎是信號處理的方面的重點只是呢,現在就先學習學習傅里葉變換吧。

   

上面這幅圖在知乎一個很著名的關於傅里葉變換的文章中的核心插圖,我覺得這幅圖很直觀的就說明了傅里葉變換的實質。時域上的東西直觀的反應到了頻域上了,很完美的結合到了一起,233333.  無數正弦波疊加,震盪的疊加的最后結果竟然是方波,同理,任何周期性函數竟然都能拆分為傅里葉級數的形式,這樣的簡介與優雅,真令人折服。

2.MATLAB對信號做頻譜分析

代碼:(1)對 f1 = Sa(2t)的頻譜分析   

 1 clear;clc;
 2 hold on;                            
 3 R=0.05;
 4 t=-1.2:R:1.2;
 5 t1 = 2*t;
 6 f1=sinc(t1);                       %Sa函數                        
 7 subplot(1,2,1),plot(t,f1)
 8 xlabel('t'),ylabel('f1')
 9 axis([-2,2,-0.3,1.2]);             %寫出Sa函數上下限
10 
11 N=1000;
12 k=-N:N;
13 W1=40;
14 W=k*W1/N;
15 F=f1*exp(-j*t'*W)*R;               %f1的傅里葉變換
16 F=real(F);                         %取F的實部
17 subplot(1,2,2),plot(W,F)
18 xlabel('W'),ylabel('F(jw)')  
View Code

結果如下圖:

     

(2)對 f2 = u(t+2) - u(t-2)的頻譜分析

 1 R=0.05;
 2 t=-3:R:3;
 3 f2=(t>=-2)-(t>=2);
 4 subplot(1,2,1),plot(t,f2)
 5 grid on;
 6 xlabel('t'),ylabel('f2')
 7 axis([-3,3,-0.5,1.5]);
 8 
 9 N=1000;k=-N:N;
10 W1=40;
11 W=k*W1/N;
12 F=f2*exp(-j*t'*W)*R;
13 F=real(F);
14 subplot(1,2,2),plot(W,F)
15 grid on;
16 xlabel('W'),ylabel('F(jw)')  
View Code

結果如下圖:

     

(3)對f3 = t[u(t+1) - u(t-1) ]的頻譜分析

 1 R=0.05;
 2 h=0.001;
 3 t=-1.2:R:1.2;
 4 y=t.*(t>=-1)-t.*(t>=1);
 5 f4=diff(y)/h;
 6 subplot(1,2,1),plot(t,y) 
 7 xlabel('t'),ylabel('y')
 8 axis([-1.2,1.2,-1.2,1.2]);
 9 
10 N=1000;
11 k=-N:N;
12 W1=40;
13 W=k*W1/N;
14 F=y*exp(-j*t'*W)*R;
15 F=real(F);
16 subplot(1,2,2),plot(W,F)
17 xlabel('W'),ylabel('F(jw)')
18 axis([-40,40,-0.06,0.06]);
View Code

結果如下圖:

     

 (4)對正弦波做FFT頻譜分析

 1 %*************************************************************************% 
 2 %                              FFT實踐及頻譜分析                           % 
 3 %*************************************************************************% 
 4 %***************正弦波****************% 
 5 fs=100;%設定采樣頻率 
 6 N=128; 
 7 n=0:N-1; 
 8 t=n/fs; 
 9 f0=10;%設定正弦信號頻率 
10 %生成正弦信號 
11 x=sin(2*pi*f0*t); 
12 figure(1); 
13 subplot(231); 
14 plot(t,x);%作正弦信號的時域波形 
15 xlabel('t'); 
16 ylabel('y'); 
17 title('正弦信號y=2*pi*10t時域波形'); 
18 grid; 
19 
20 %進行FFT變換並做頻譜圖 
21 y=fft(x,N);%進行fft變換 
22 mag=abs(y);%求幅值 
23 f=(0:length(y)-1)'*fs/length(y);%進行對應的頻率轉換 
24 figure(1); 
25 subplot(232); 
26 plot(f,mag);%做頻譜圖 
27 axis([0,100,0,80]); 
28 xlabel('頻率(Hz)'); 
29 ylabel('幅值'); 
30 title('正弦信號y=2*pi*10t幅頻譜圖N=128'); 
31 grid; 
32 
33 %求均方根譜 
34 sq=abs(y); 
35 figure(1); 
36 subplot(233); 
37 plot(f,sq); 
38 xlabel('頻率(Hz)'); 
39 ylabel('均方根譜'); 
40 title('正弦信號y=2*pi*10t均方根譜'); 
41 grid; 
42 
43 %求功率譜 
44 power=sq.^2; 
45 figure(1); 
46 subplot(234); 
47 plot(f,power); 
48 xlabel('頻率(Hz)'); 
49 ylabel('功率譜'); 
50 title('正弦信號y=2*pi*10t功率譜'); 
51 grid; 
52 
53 %求對數譜 
54 ln=log(sq); 
55 figure(1); 
56 subplot(235); 
57 plot(f,ln); 
58 xlabel('頻率(Hz)'); 
59 ylabel('對數譜'); 
60 title('正弦信號y=2*pi*10t對數譜'); 
61 grid; 
62 
63 %用IFFT恢復原始信號 
64 xifft=ifft(y); 
65 magx=real(xifft); 
66 ti=[0:length(xifft)-1]/fs; 
67 figure(1); 
68 subplot(236); 
69 plot(ti,magx); 
70 xlabel('t'); 
71 ylabel('y'); 
72 title('通過IFFT轉換的正弦信號波形'); 
73 grid; 
View Code

執行結果如下圖:

    

 

(5)對矩形波做FFT頻譜分析

 1 %****************2.矩形波****************% 
 2 fs=10;%設定采樣頻率 
 3 t=-5:0.1:5; 
 4 x=rectpuls(t,2); 
 5 x=x(1:99); 
 6 figure(1); 
 7 subplot(231); plot(t(1:99),x);%作矩形波的時域波形 
 8 xlabel('t'); 
 9 ylabel('y'); 
10 title('矩形波時域波形'); 
11 grid; 
12 
13 %進行FFT變換並做頻譜圖 
14 y=fft(x);%進行fft變換 
15 mag=abs(y);%求幅值 
16 f=(0:length(y)-1)'*fs/length(y);%進行對應的頻率轉換 
17 figure(1); 
18 subplot(232); 
19 plot(f,mag);%做頻譜圖 
20 xlabel('頻率(Hz)'); 
21 ylabel('幅值'); 
22 title('矩形波幅頻譜圖'); 
23 grid; 
24 
25 %求均方根譜 
26 sq=abs(y); 
27 figure(1); 
28 subplot(233); 
29 plot(f,sq); 
30 xlabel('頻率(Hz)'); 
31 ylabel('均方根譜'); 
32 title('矩形波均方根譜'); 
33 grid; 
34 
35 %求功率譜 
36 power=sq.^2; 
37 figure(1); 
38 subplot(234); 
39 plot(f,power); 
40 xlabel('頻率(Hz)'); 
41 ylabel('功率譜'); 
42 title('矩形波功率譜'); 
43 grid; 
44 
45 %求對數譜 
46 ln=log(sq); 
47 figure(1); 
48 subplot(235); 
49 plot(f,ln); 
50 xlabel('頻率(Hz)'); 
51 ylabel('對數譜'); 
52 title('矩形波對數譜'); 
53 grid; 
54 
55 %用IFFT恢復原始信號 
56 xifft=ifft(y); 
57 magx=real(xifft); 
58 ti=[0:length(xifft)-1]/fs; 
59 figure(1); 
60 subplot(236); 
61 plot(ti,magx); 
62 xlabel('t'); 
63 ylabel('y'); 
64 title('通過IFFT轉換的矩形波波形'); 
65 grid; 
View Code

執行結果如下圖:

       

(6)對白噪聲做頻譜分析

 1 %****************3.白噪聲****************% 
 2 fs=10;%設定采樣頻率 
 3 t=-5:0.1:5; 
 4 x=zeros(1,100); 
 5 x(50)=100000; 
 6 figure(1); 
 7 subplot(231); 
 8 plot(t(1:100),x);%作白噪聲的時域波形 
 9 xlabel('t'); 
10 ylabel('y'); 
11 title('白噪聲時域波形'); 
12 grid; 
13 
14 %進行FFT變換並做頻譜圖 
15 y=fft(x);  %進行fft變換 
16 mag=abs(y);%求幅值 
17 f=(0:length(y)-1)'*fs/length(y);%進行對應的頻率轉換 
18 figure(1); 
19 subplot(232); 
20 plot(f,mag);%做頻譜圖 
21 xlabel('頻率(Hz)'); 
22 ylabel('幅值'); 
23 title('白噪聲幅頻譜圖'); 
24 grid; 
25 
26 %求均方根譜 
27 sq=abs(y); 
28 figure(1); 
29 subplot(233); 
30 plot(f,sq); 
31 xlabel('頻率(Hz)'); 
32 ylabel('均方根譜'); 
33 title('白噪聲均方根譜'); 
34 grid; 
35 
36 %求功率譜 
37 power=sq.^2; 
38 figure(1); 
39 subplot(234); 
40 plot(f,power); 
41 xlabel('頻率(Hz)'); 
42 ylabel('功率譜'); 
43 title('白噪聲功率譜'); 
44 grid; 
45 
46 %求對數譜 
47 ln=log(sq); 
48 figure(1); 
49 subplot(235); 
50 plot(f,ln); 
51 xlabel('頻率(Hz)'); 
52 ylabel('對數譜'); 
53 title('白噪聲對數譜'); 
54 grid; 
55 
56 %用IFFT恢復原始信號 
57 xifft=ifft(y); 
58 magx=real(xifft); 
59 ti=[0:length(xifft)-1]/fs; 
60 figure(1); 
61 subplot(236); 
62 plot(ti,magx); 
63 xlabel('t'); 
64 ylabel('y'); 
65 title('通過IFFT轉換的白噪聲波形'); 
66 grid;   
View Code

執行結果如下:

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM