數(shù)字信號MATLAB仿真實驗

數(shù)字信號MATLAB仿真實驗實驗?zāi)康耐ㄟ^本次實驗熟悉數(shù)字信號處理的方波和三角波的產(chǎn)生對FFT加強理解和直觀地認識理解并掌握頻譜分析和諧波分析對噪聲對信號的影響認識實驗內(nèi)容自編程序得到一個方波信號（f=50Hz,幅值為1，，各半個周期），對其一個周期分別采樣256點和1024點，利用光盤所附的基于Matlab語言所編FFT程序或自編FFT程序做諧波分析分析，并與理論分析結(jié)果對照（注意FFT計算頻譜與諧波分析的區(qū)別）。對三角波信號（可以由方波信號求導(dǎo)得到）重復(fù)作業(yè)一的各項要求。對一、二信號疊加一個白噪聲信號（均值為零，方差為0.2）所構(gòu)成的隨機信號用FFT進行頻譜分析對以上結(jié)果進行討論實驗程序方波分析cclearT=0.02,l=0.5,hu=1,hd=-1;%T是周期，l為占空比，hu為上幅度，hd為下幅度t=[0:0.001:2];%顯示連續(xù)圖形的采樣點fori=1:length(t)%生成方波ifrem(t(i),T)<T*ly(i)=hu;elsey(i)=hd;endendn=0:0.02/256:(0.02-0.02/256);%256采樣點fori=1:length(n)%生成方波取256采樣點值ifrem(n(i),T)<T*lyn(i)=hu;elseyn(i)=hd;endendn2=0:0.02/1024:(0.02-0.02/1024);fori=1:length(n2)%生成方波取1024采樣點值ifrem(n2(i),T)<T*lyn2(i)=hu;elseyn2(i)=hd;endendfigure(1)subplot(3,1,1)stairs(t,y); %畫方波xlabel('時間t'),ylabel('幅度y'),gridon;subplot(3,1,2); %畫256采樣圖stem(n,yn,'.');xlabel('采樣時間t'),ylabel('幅度y'),gridon;subplot(3,1,3); %畫1024采樣圖stem(n2,yn2,'.');xlabel('采樣時間t'),ylabel('幅度y'),gridon;ynfft=fft(yn); %256點FFTynfft2=fft(yn2); %1024點FFTmagyn=abs(ynfft); %256點幅度譜，對數(shù)幅度譜，相位譜lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2); %1024點幅度譜，對數(shù)幅度譜，相位譜lgmagyn2=20*log10(abs(ynfft2));argyn2=angle(ynfft2);nf=0:2*pi/256:(2*pi-2*pi/256);%計算FFT各點頻率nf2=0:2*pi/1024:(2*pi-2*pi/1024);ind1=find(magyn);%取fft后的幅值非零頻率分量magynxb=magyn(:,ind1);ind2=find(magyn2);magynxb2=magyn2(:,ind2);ifind1(1)~=1%判斷是否有直流分量，并求基波頻率 %調(diào)整各各次諧波的存儲在相應(yīng)位置bf1=gcd(ind1(1)-1,ind1(2)-1)*2*pi/256;elsebf1=gcd(ind1(2)-1,ind1(3)-1)*2*pi/256;magynxb=magyn(:,2:ind1);ind11=ind1(2:length(ind1))-1;endifind2(1)~=1bf2=gcd(ind2(1)-1,ind2(2)-1)*2*pi/1024;elsebf2=gcd(ind2(2)-1,ind2(3)-1)*2*pi/256;magynxb2=magyn2(2:1024);%去掉直流分量ind2=ind2(2:length(ind2))-1endfigure(2) %畫各諧波的幅度譜subplot(2,1,1)stem(ind1-1,magynxb,'.');xlabel('諧波頻率'),ylabel('幅值'),gridon;subplot(2,1,2)stem(ind2-1,magynxb2,'.');xlabel('諧波頻率'),ylabel('幅值'),gridon;learT=0.02,l=0.5,hu=1,hd=-1;%T是周期，l為占空比，hu為上幅度，hd為下幅度t=[0:0.001:2];%顯示連續(xù)圖形的采樣點fori=1:length(t)%生成方波ifrem(t(i),T)<T*ly(i)=hu;elsey(i)=hd;endendn=0:0.02/256:(0.02-0.02/256);%256采樣點fori=1:length(n)%生成方波取256采樣點值ifrem(n(i),T)<T*lyn(i)=hu;elseyn(i)=hd;endendn2=0:0.02/1024:(0.02-0.02/1024);fori=1:length(n2)%生成方波取1024采樣點值ifrem(n2(i),T)<T*lyn2(i)=hu;elseyn2(i)=hd;endendfigure(1)subplot(3,1,1)stairs(t,y);%畫方波xlabel('時間t'),ylabel('幅度y'),gridon;subplot(3,1,2);%畫256采樣圖stem(n,yn,'.');xlabel('采樣時間t'),ylabel('幅度y'),gridon;subplot(3,1,3);%畫1024采樣圖stem(n2,yn2,'.');xlabel('采樣時間t'),ylabel('幅度y'),gridon;ynfft=fft(yn);%256點FFTynfft2=fft(yn2);%1024點FFTmagyn=abs(ynfft);%256點幅度譜，對數(shù)幅度譜，相位譜lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2);%1024點幅度譜，對數(shù)幅度譜，相位譜lgmagyn2=20*log10(abs(ynfft2));argyn2=angle(ynfft2);nf=0:2*pi/256:(2*pi-2*pi/256);%計算FFT各點頻率nf2=0:2*pi/1024:(2*pi-2*pi/1024);ind1=find(magyn);%取fft后的幅值非零頻率分量magynxb=magyn(:,ind1);ind2=find(magyn2);magynxb2=magyn2(:,ind2);ifind1(1)~=1%判斷是否有直流分量，并求基波頻率bf1=gcd(ind1(1)-1,ind1(2)-1)*2*pi/256;elseind1bf1=gcd(ind1(2)-1,ind1(3)-1)*2*pi/256;magynxb=magyn(:,2:ind1);endifind2(1)~=1bf2=gcd(ind2(1)-1,ind2(2)-1)*2*pi/1024;elsebf2=gcd(ind2(2)-1,ind2(3)-1)*2*pi/256;magynxb2=magyn2(2:1024);%去掉直流分量endfigure(2)%畫各諧波的幅度譜subplot(2,1,1)stem(ind1,magynxb,'.');xlabel('頻率'),ylabel('幅值'),gridon;subplot(2,1,2)stem(ind2,magynxb2,'.');xlabel('頻率'),ylabel('幅值'),gridon;三角波分析clearT=0.02,l=0.5,hu=1,hd=-1;%T是周期，l為占空比，hu為上幅度，hd為下幅度t=[0:0.0001:2];%顯示連續(xù)圖形的采樣點fori=1:length(t)%生成三角波ifrem(t(i),T)<T*ly(i)=-(hu-hd)/(T*l)*rem(t(i),T)+hu;elsey(i)=(hu-hd)/(T*(1-l))*(rem(t(i),T)-l*T)+hd;endendn=0:0.02/256:(0.02-0.02/256);%256采樣點fori=1:length(n)%生成256點采樣值ifrem(n(i),T)<T*lyn(i)=-(hu-hd)/(T*l)*rem(n(i),T)+hu;elseyn(i)=(hu-hd)/(T*(1-l))*(rem(n(i),T)-l*T)+hd;endendn2=0:0.02/1024:(0.02-0.02/1024);fori=1:length(n2)%生成1024點采樣值ifrem(n2(i),T)<T*lyn2(i)=-(hu-hd)/(T*l)*rem(n2(i),T)+hu;elseyn2(i)=(hu-hd)/(T*(1-l))*(rem(n2(i),T)-l*T)+hd;endendfigure(1)subplot(3,1,1)plot(t,y);%畫三角波xlabel('時間t'),ylabel('幅度y'),gridon;subplot(3,1,2);%畫256采樣圖stem(n,yn,'.');xlabel('采樣時間t'),ylabel('幅度y'),gridon;subplot(3,1,3);%畫1024采樣圖stem(n2,yn2,'.');xlabel('采樣時間t'),ylabel('幅度y'),gridon;ynfft=fft(yn);%256點FFTynfft2=fft(yn2);%1024點FFTmagyn=abs(ynfft);%256點幅度譜，對數(shù)幅度譜，相位譜lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2);%1024點幅度譜，對數(shù)幅度譜，相位譜lgmagyn2=20*log10(abs(ynfft2));argyn2=angle(ynfft2);nf=0:2*pi/256:(2*pi-2*pi/256);%計算FFT各點頻率nf2=0:2*pi/1024:(2*pi-2*pi/1024);ind1=find(magyn);%取fft后的幅值非零頻率分量magynxb=magyn(:,ind1);ind2=find(magyn2);magynxb2=magyn2(:,ind2);ifind1(1)~=1%判斷是否有直流分量，并求基波頻率bf1=gcd(ind1(1)-1,ind1(2)-1)*2*pi/256;elsebf1=gcd(ind1(2)-1,ind1(3)-1)*2*pi/256;magynxb=magyn(2:length(magyn));ind11=ind1(2:length(ind1));endifind2(1)~=1bf2=gcd(ind2(1)-1,ind2(2)-1)*2*pi/1024;elsebf2=gcd(ind2(2)-1,ind2(3)-1)*2*pi/256;magynxb2=magyn2(2:length(magyn2));%去掉直流分量ind22=ind2(2:length(ind2));endfigure(2)%畫各諧波的幅度譜subplot(2,1,1)stem(ind11-1,magynxb,'.');xlabel('頻率'),ylabel('幅值'),gridon;subplot(2,1,2)stem(ind22-1,magynxb2,'.');xlabel('頻率'),ylabel('幅值'),gridon;疊加噪聲分析疊加了噪聲的程序TclearT=0.02,l=0.5,hu=1,hd=-1;%T是周期，l為占空比，hu為上幅度，hd為下幅度t=[0:0.0001:2];%顯示連續(xù)圖形的采樣點fori=1:length(t)%生成方波ifrem(t(i),T)<T*lyrect(i)=hu;elseyrect(i)=hd;endendfori=1:length(t)%生成三角波ifrem(t(i),T)<T*lytri(i)=-(hu-hd)/(T*l)*rem(t(i),T)+hu;elseytri(i)=(hu-hd)/(T*(1-l))*(rem(t(i),T)-l*T)+hd;endendnoise=sqrt(0.2)*randn(1,length(t));%noise=sqrt(12*0.2)*rand(1,length(t))y=yrect+noise;y2=ytri+noise;n=0:0.02/1024:(0.02-0.02/1024);%1024采樣點n2=nfori=1:length(n)%生成方波取1024采樣點值ifrem(n(i),T)<T*lyn(i)=hu;elseyn(i)=hd;endendfori=1:length(n2)%生成1024點三角波采樣值ifrem(n2(i),T)<T*lyn2(i)=-(hu-hd)/(T*l)*rem(n2(i),T)+hu;elseyn2(i)=(hu-hd)/(T*(1-l))*(rem(n2(i),T)-l*T)+hd;endendn=[0:0.02/1024:0.02-0.02/1024];n2=n;noisen=sqrt(0.2)*randn(1,length(n));%noise=sqrt(12*0.2)*rand(1,length(n))yn=yn+noisen;yn2=yn2+noisen;figure(1)subplot(4,2,1)stairs(t,yrect)ylabel('方波圖'),gridon;subplot(4,2,5)plot(t,y);ylabel('疊加了噪聲的方波圖'),gridon;subplot(4,2,7);stem(n,yn,'.')ylabel('疊加了噪聲的方波圖的抽樣'),gridon;subplot(4,2,3)plot(t,noise)ylabel('正態(tài)分布噪聲'),gridon;subplot(4,2,4)ylabel('噪聲的分布檢測'),gridon;hist(noise,50)subplot(4,2,2)stairs(t,ytri)ylabel('三角波圖'),gridon;subplot(4,2,6)plot(t,y2)ylabel('疊加了噪聲的三角波圖'),gridon;subplot(4,2,8)ylabel('疊加了噪聲的三角波圖的抽樣'),gridon;stem(n2,yn2,'.')ynfft=fft(yn