利用Matlab繪制正弦信號(hào)的頻譜圖并做相關(guān)分析

1、利用Matlab繪制正弦信號(hào)的頻譜圖并做相關(guān)分析一、作業(yè)要求：1、信號(hào)可變（信號(hào)的賦值、相位、頻率可變）；2、采樣頻率fs可變；3、加各種不同的窗函數(shù)并分析其影響；4、頻譜校正；5、頻譜細(xì)化。二、采用matlab編寫(xiě)如下程序：clear;clf;fs=100;N=1024; %采樣頻率和數(shù)據(jù)點(diǎn)數(shù)A=20;B=30;C=0.38;n=0:N-1;t=n/fs; %時(shí)間序列x=A*sin(2*pi*B*t+C); %信號(hào)y=fft(x,N); %對(duì)信號(hào)進(jì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;%兩種信號(hào)疊加，x=A*sin(2*pi*B*t+C)+2*A*sin(2*pi*1.5*B*t+2.5*C); %信號(hào)y=fft(x,N); %對(duì)信號(hào)進(jì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，兩種信號(hào)疊加');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;%改變采樣點(diǎn)數(shù)N=128N=128;n=0:N-1;t=n/fs; %時(shí)間序列x=A*sin(2*pi*B*t+C); %信號(hào)y=fft(x,N); %對(duì)信號(hào)進(jì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時(shí)的頻譜fs=400;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信號(hào)y=fft(x,N); %對(duì)信號(hào)進(jì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ù)點(diǎn)數(shù)n=0:N-1;t=n/fs; %時(shí)間序列x=A*sin(2*pi*B*t+C); %信號(hào)window=triang(N);%生成三角窗函數(shù)x=x.*window'%加窗函數(shù)y=fft(x,N); %對(duì)信號(hào)進(jì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); %信號(hào)window=hamming(N);%生成海明窗函數(shù)x=x.*window'%加窗函數(shù)y=fft(x,N); %對(duì)信號(hào)進(jì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); %信號(hào)window=hanning(N);%生成漢寧窗函數(shù)x=x.*window'%加窗函數(shù)y=fft(x,N); %對(duì)信號(hào)進(jì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;三、運(yùn)行結(jié)果如下：四、分析與結(jié)論：1）從所做圖像可以看出，信號(hào)的幅值均小于真實(shí)值，說(shuō)明在截?cái)嘈盘?hào)時(shí)存在泄露。2）從圖1和圖圖2取相同的采樣頻率fs=100和數(shù)據(jù)點(diǎn)數(shù)N=1024，不同的是圖2采用兩種不同賦值和頻率的正弦信號(hào)疊加，從圖中可以看出，圖2可以明顯的看出含有兩種不同的頻率成分的信號(hào)，幅值也不相同，由此可以看出，不同頻率的正弦信號(hào)疊加，在頻域當(dāng)中互相分離，互不影響

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

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

12、ning(N/2)');%對(duì)信號(hào)做N/2點(diǎn)FFT變換Y1=abs(y1(1:N/2)/N*2);%第一段信號(hào)幅值Y2=abs(y2(1:N/4)/N*4);%第二段信號(hào)幅值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調(diào)整到（-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