數(shù)字信號(hào)MATLAB仿真實(shí)驗(yàn)

數(shù)字信號(hào)MATLAB仿真實(shí)驗(yàn)實(shí)驗(yàn)?zāi)康耐ㄟ^本次實(shí)驗(yàn)熟悉數(shù)字信號(hào)處理的方波和三角波的產(chǎn)生對(duì)FFT加強(qiáng)理解和直觀地認(rèn)識(shí)理解并掌握頻譜分析和諧波分析對(duì)噪聲對(duì)信號(hào)的影響認(rèn)識(shí)實(shí)驗(yàn)內(nèi)容自編程序得到一個(gè)方波信號(hào)（f=50Hz,幅值為1，，各半個(gè)周期），對(duì)其一個(gè)周期分別采樣256點(diǎn)和1024點(diǎn)，利用光盤所附的基于Matlab語言所編FFT程序或自編FFT程序做諧波分析分析，并與理論分析結(jié)果對(duì)照（注意FFT計(jì)算頻譜與諧波分析的區(qū)別）。對(duì)三角波信號(hào)（可以由方波信號(hào)求導(dǎo)得到）重復(fù)作業(yè)一的各項(xiàng)要求。對(duì)一、二信號(hào)疊加一個(gè)白噪聲信號(hào)（均值為零，方差為0.2）所構(gòu)成的隨機(jī)信號(hào)用FFT進(jìn)行頻譜分析對(duì)以上結(jié)果進(jìn)行討論實(shí)驗(yàn)程序方波分析cclearT=0.02,l=0.5,hu=1,hd=-1;%T是周期，l為占空比，hu為上幅度，hd為下幅度t=[0:0.001:2];%顯示連續(xù)圖形的采樣點(diǎn)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采樣點(diǎn)fori=1:length(n)%生成方波取256采樣點(diǎn)值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采樣點(diǎn)值ifrem(n2(i),T)<T*lyn2(i)=hu;elseyn2(i)=hd;endendfigure(1)subplot(3,1,1)stairs(t,y); %畫方波xlabel('時(shí)間t'),ylabel('幅度y'),gridon;subplot(3,1,2); %畫256采樣圖stem(n,yn,'.');xlabel('采樣時(shí)間t'),ylabel('幅度y'),gridon;subplot(3,1,3); %畫1024采樣圖stem(n2,yn2,'.');xlabel('采樣時(shí)間t'),ylabel('幅度y'),gridon;ynfft=fft(yn); %256點(diǎn)FFTynfft2=fft(yn2); %1024點(diǎn)FFTmagyn=abs(ynfft); %256點(diǎn)幅度譜，對(duì)數(shù)幅度譜，相位譜lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2); %1024點(diǎn)幅度譜，對(duì)數(shù)幅度譜，相位譜lgmagyn2=20*log10(abs(ynfft2));argyn2=angle(ynfft2);nf=0:2*pi/256:(2*pi-2*pi/256);%計(jì)算FFT各點(diǎn)頻率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)整各各次諧波的存儲(chǔ)在相應(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ù)圖形的采樣點(diǎn)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采樣點(diǎn)fori=1:length(n)%生成方波取256采樣點(diǎn)值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采樣點(diǎn)值ifrem(n2(i),T)<T*lyn2(i)=hu;elseyn2(i)=hd;endendfigure(1)subplot(3,1,1)stairs(t,y);%畫方波xlabel('時(shí)間t'),ylabel('幅度y'),gridon;subplot(3,1,2);%畫256采樣圖stem(n,yn,'.');xlabel('采樣時(shí)間t'),ylabel('幅度y'),gridon;subplot(3,1,3);%畫1024采樣圖stem(n2,yn2,'.');xlabel('采樣時(shí)間t'),ylabel('幅度y'),gridon;ynfft=fft(yn);%256點(diǎn)FFTynfft2=fft(yn2);%1024點(diǎn)FFTmagyn=abs(ynfft);%256點(diǎn)幅度譜，對(duì)數(shù)幅度譜，相位譜lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2);%1024點(diǎn)幅度譜，對(duì)數(shù)幅度譜，相位譜lgmagyn2=20*log10(abs(ynfft2));argyn2=angle(ynfft2);nf=0:2*pi/256:(2*pi-2*pi/256);%計(jì)算FFT各點(diǎn)頻率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ù)圖形的采樣點(diǎn)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采樣點(diǎn)fori=1:length(n)%生成256點(diǎn)采樣值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點(diǎn)采樣值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('時(shí)間t'),ylabel('幅度y'),gridon;subplot(3,1,2);%畫256采樣圖stem(n,yn,'.');xlabel('采樣時(shí)間t'),ylabel('幅度y'),gridon;subplot(3,1,3);%畫1024采樣圖stem(n2,yn2,'.');xlabel('采樣時(shí)間t'),ylabel('幅度y'),gridon;ynfft=fft(yn);%256點(diǎn)FFTynfft2=fft(yn2);%1024點(diǎn)FFTmagyn=abs(ynfft);%256點(diǎn)幅度譜，對(duì)數(shù)幅度譜，相位譜lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2);%1024點(diǎn)幅度譜，對(duì)數(shù)幅度譜，相位譜lgmagyn2=20*log10(abs(ynfft2));argyn2=angle(ynfft2);nf=0:2*pi/256:(2*pi-2*pi/256);%計(jì)算FFT各點(diǎn)頻率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ù)圖形的采樣點(diǎn)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采樣點(diǎn)n2=nfori=1:length(n)%生成方波取1024采樣點(diǎn)值ifrem(n(i),T)<T*lyn(i)=hu;elseyn(i)=hd;endendfori=1:length(n2)%生成1024點(diǎn)三角波采樣值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('噪聲的分布檢測(cè)'),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