利用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 編寫如下程序：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

2、,1),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(頻率/itHz);ylabel(振幅);title(圖 2：fs=100,

3、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 混入噪聲);無grid 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)行傅里葉變

4、換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;%改變采樣頻率為 200Hz 時(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

5、=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;%采樣頻率和數(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;%幅值處

6、理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：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=hamming(N);%生成海明窗函數(shù)x=x.*window;%加窗函數(shù)無y=fft(x,N);%對(duì)信號(hào)進(jìn)行快速傅里葉變換yy=abs(y);%求取傅里葉變換的振幅yy=yy*2

7、/N;%幅值處理f=n*fs/N;subplot(3,3,7),plot(f(1:N/2.56),1.852*yy(1:N/2.56); %繪出隨頻率變化的振幅xlabel(頻率/itHz);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);%求取傅里葉變換的振幅y

8、y=yy*2/N;%幅值處理f=n*fs/N;subplot(3,3,8),plot(f(1:N/2.56),2*yy(1:N/2.56); %繪出隨頻率變化的振幅xlabel(頻率/itHz);ylabel(振幅);title(圖 8：fs=100，N=1024,加漢寧窗函數(shù));grid on;三、運(yùn)行結(jié)果如下：無四、分析與結(jié)論：1）從所做圖像可以看出，信號(hào)的幅值均小于真實(shí)值，說明在截?cái)嘈盘?hào)時(shí)存在泄露。2）從圖 1 和圖圖 2 取相同的采樣頻率 fs=100 和數(shù)據(jù)點(diǎn)數(shù) N=1024，不同的是圖 2 采用兩種不同賦值和頻率的正弦信號(hào)疊加，從圖中可以看出，圖 2 可以明顯的看出含有兩種不同的頻

9、率成分的信號(hào)，幅值也不相同，由此可以看出，不同頻率的正弦信號(hào)疊加，在頻域當(dāng)中互相分離，互不影響。3）從圖 1 和圖圖 2 可以看出，整個(gè)頻譜圖是以 fs/2 頻率為對(duì)稱軸的。由此可以知道傅里葉變換數(shù)據(jù)的對(duì)稱性。因此在用傅里葉變換做頻譜分析時(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ù)有

10、關(guān)。一定范圍內(nèi)采樣點(diǎn)數(shù)越多，信號(hào)的幅值越接近真實(shí)值。6）圖 5 改變采樣頻率觀察不同采樣頻率對(duì)信號(hào)的影響，當(dāng)采樣頻率太小時(shí)，譜線的尾部發(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

11、*B*t+C); %正弦信號(hào)y1=fft(x.*hanning(N);%對(duì)信號(hào)做 N 點(diǎn) FFT 變換y2=fft(x(1:N/2).*hanning(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=ang

12、le(y1(k1);phase2=angle(y2(k2);Ano=Y1Amax*2fno=(k1-1)*fs/N%未校正頻率phaseno=phase1*180/pi%未校正相角無delt=mod(phase1-phase2,2*pi);%將 delt 調(diào)整到（-pi,pi）之間if deltpidelt1=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