利用Matlab繪制正弦信號的頻譜圖并做相關分析

1、利用Matlab繪制正弦信號的頻譜圖并做相關分析一、作業(yè)要求：1、信號可變（信號的賦值、相位、頻率可變）；2、采樣頻率fs可變；3、加各種不同的窗函數(shù)并分析其影響；4、頻譜校正；5、頻譜細化。二、采用matlab編寫如下程序：clear;clf;fs=100;N=1024; %采樣頻率和數(shù)據(jù)點數(shù)A=20;B=30;C=0.38;n=0:N-1;t=n/fs; %時間序列x=A*sin(2*pi*B*t+C); %信號y=fft(x,N); %對信號進行傅里葉變換yy=abs(y); %求得傅里葉變換后的振幅yy=yy*2/N; %幅值處理f=n*fs/N; %頻率序列subplot(3,3,1

2、),plot(f,yy); %繪出隨頻率變化的振幅xlabel('頻率/itHz');ylabel('振幅');title('圖1：fs=100，N=1024');grid on;%兩種信號疊加，x=A*sin(2*pi*B*t+C)+2*A*sin(2*pi*1.5*B*t+2.5*C); %信號y=fft(x,N); %對信號進行傅里葉變換yy=abs(y); %求得傅里葉變換后的振幅yy=yy*2/N; %幅值處理f=n*fs/N; %頻率序列subplot(3,3,2),plot(f,yy); %繪出隨頻率變化的振幅xlabel('

3、;頻率/itHz');ylabel('振幅');title('圖2：fs=100,N=1024，兩種信號疊加');grid on;%加噪聲之后的圖像x=A*sin(2*pi*B*t+C)+28*randn(size(t);y=fft(x,N);yy=abs(y);yy=yy*2/N; %幅值處理subplot(3,3,3),plot(f(1:N/2.56),yy(1:N/2.56);xlabel('頻率/itHz');ylabel('振幅');title('圖3：fs=100,N=1024混入噪聲');gr

4、id on;%改變采樣點數(shù)N=128N=128;n=0:N-1;t=n/fs; %時間序列x=A*sin(2*pi*B*t+C); %信號y=fft(x,N); %對信號進行傅里葉變換yy=abs(y); %求得傅里葉變換后的振幅yy=yy*2/N; %幅值處理f=n*fs/N; %頻率序列subplot(3,3,4),plot(f(1:N/2.56),yy(1:N/2.56); %繪出隨頻率變化的振幅xlabel('頻率/itHz');ylabel('振幅');title('圖4：fs=100，N=128');grid on;%改變采樣頻率為2

5、00Hz時的頻譜fs=400;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信號y=fft(x,N); %對信號進行快速傅里葉變換yy=abs(y); %求取傅里葉變換的振幅yy=yy*2/N; %幅值處理f=n*fs/N;subplot(3,3,5),plot(f(1:N/2.56),yy(1:N/2.56); %繪出隨頻率變化的振幅xlabel('頻率/itHz');ylabel('振幅');title('圖5：fs=400，N=1024');grid on;%加三角窗函數(shù)fs=100;N=1024;

6、 %采樣頻率和數(shù)據(jù)點數(shù)n=0:N-1;t=n/fs; %時間序列x=A*sin(2*pi*B*t+C); %信號window=triang(N);%生成三角窗函數(shù)x=x.*window'%加窗函數(shù)y=fft(x,N); %對信號進行傅里葉變換yy=abs(y); %求得傅里葉變換后的振幅yy=yy*2/N; %幅值處理f=n*fs/N; %頻率序列subplot(3,3,6),plot(f(1:N/2.56),2*yy(1:N/2.56); %繪出隨頻率變化的振幅xlabel('頻率/itHz');ylabel('振幅');title('圖6：f

7、s=100,N=1024,加三角窗函數(shù)');grid on;%加海明窗函數(shù)后的頻譜fs=100;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信號window=hamming(N);%生成海明窗函數(shù)x=x.*window'%加窗函數(shù)y=fft(x,N); %對信號進行快速傅里葉變換yy=abs(y); %求取傅里葉變換的振幅yy=yy*2/N; %幅值處理f=n*fs/N;subplot(3,3,7),plot(f(1:N/2.56),1.852*yy(1:N/2.56); %繪出隨頻率變化的振幅xlabel('頻率/itHz&

8、#39;);ylabel('振幅');title('圖7：fs=100，N=1024,加海明窗函數(shù)');grid on;%加漢寧窗函數(shù)后的頻譜fs=100;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信號window=hanning(N);%生成漢寧窗函數(shù)x=x.*window'%加窗函數(shù)y=fft(x,N); %對信號進行快速傅里葉變換yy=abs(y); %求取傅里葉變換的振幅yy=yy*2/N; %幅值處理f=n*fs/N;subplot(3,3,8),plot(f(1:N/2.56),2*yy(1:N/

9、2.56); %繪出隨頻率變化的振幅xlabel('頻率/itHz');ylabel('振幅');title('圖8：fs=100，N=1024,加漢寧窗函數(shù)');grid on;三、運行結果如下：四、分析與結論：1）從所做圖像可以看出，信號的幅值均小于真實值，說明在截斷信號時存在泄露。2）從圖1和圖圖2取相同的采樣頻率fs=100和數(shù)據(jù)點數(shù)N=1024，不同的是圖2采用兩種不同賦值和頻率的正弦信號疊加，從圖中可以看出，圖2可以明顯的看出含有兩種不同的頻率成分的信號，幅值也不相同，由此可以看出，不同頻率的正弦信號疊加，在頻域當中互相分離，互不影響

10、。3）從圖1和圖圖2可以看出，整個頻譜圖是以fs/2頻率為對稱軸的。由此可以知道傅里葉變換數(shù)據(jù)的對稱性。因此在用傅里葉變換做頻譜分析時，我們只需做出前一半頻譜圖即可。4）圖3為混入噪聲之后的頻譜，可以看出噪聲分布在整個頻率軸上，并且由于噪聲中含有與原信號頻率相同的成分，疊加之后導致幅值增加。加大噪聲的幅值之后，將分辨不出原信號的頻率。5）圖4減少了數(shù)據(jù)點數(shù)，N=128，與圖1相比較，采用128點和1024點的相同頻率的振幅是有不同的表現(xiàn)值。因此振幅的大小與所用采樣點數(shù)有關。一定范圍內采樣點數(shù)越多，信號的幅值越接近真實值。6）圖5改變采樣頻率觀察不同采樣頻率對信號的影響，當采樣頻率太小時，譜線的

11、尾部發(fā)生混疊現(xiàn)象，當采樣頻率太大時，頻率的分辨率較低，不利于采樣，根據(jù)采樣定理，采樣頻率必須大于信號最高頻率的2倍，通常采用35倍。7）圖6、7、8分別對信號添加了三角窗、海明窗和漢寧窗，圖1比較，加窗之后信號的幅值更加接近真實值。而且使得圖像的旁瓣減小，信號的能量相對集中。五、采用相位差法進行頻譜校正校正程序代碼：clear;clf;fs=100;N=1024;n=0:N-1;t=n/fs;A=20;B=30;C=0.38;x=A*sin(2*pi*B*t+C); %正弦信號y1=fft(x.*hanning(N)');%對信號做N點FFT變換y2=fft(x(1:N/2).*han

12、ning(N/2)');%對信號做N/2點FFT變換Y1=abs(y1(1:N/2)/N*2);%第一段信號幅值Y2=abs(y2(1:N/4)/N*4);%第二段信號幅值f=(1:N/2)*fs/N;subplot(2,1,1);plot(f,2*Y1);xlabel('頻率/itHz');ylabel('振幅/A');title('加漢寧窗校正前');grid on;Y1Amax,k1=max(Y1);Y2Amax,k2=max(Y2);phase1=angle(y1(k1);phase2=angle(y2(k2);Ano=Y1Ama

13、x*2fno=(k1-1)*fs/N %未校正頻率phaseno=phase1*180/pi %未校正相角delt=mod(phase1-phase2,2*pi);%將delt調整到（-pi,pi）之間if delt<-pi delt1=delt+2*pi;elseif delt>pi delt1=delt-2*pi;else delt1=delt;enddeltf=2*(k2-1)-(k1-1)-2*delt1/pi;Yyes=zeros(1,N/2);Ayes=2/sinc(deltf)*Y1Amax*(1-deltf2)Yyes(k2)=Ayes;fyes=(k1-1-deltf)*fs/N %校正后頻率phaseyes=(phase1+deltf*pi)*180/pi %校正后相位f(k2