数字信号MATLAB仿真实验

数字信号MATLAB仿真实验实验目的通过本次实验熟悉数字信号处理的方波和三角波的产生对FFT加强理解和直观地认识理解并掌握频谱分析和谐波分析对噪声对信号的影响认识实验内容自编程序得到一个方波信号（f=50Hz,幅值为1，，各半个周期），对其一个周期分别采样256点和1024点，利用光盘所附的基于Matlab语言所编FFT程序或自编FFT程序做谐波分析分析，并与理论分析结果对照（注意FFT计算频谱与谐波分析的区别）。对三角波信号（可以由方波信号求导得到）重复作业一的各项要求。对一、二信号叠加一个白噪声信号（均值为零，方差为0.2）所构成的随机信号用FFT进行频谱分析对以上结果进行讨论实验程序方波分析cclearT=0.02,l=0.5,hu=1,hd=-1;%T是周期，l为占空比，hu为上幅度，hd为下幅度t=[0:0.001:2];%显示连续图形的采样点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点幅度谱，对数幅度谱，相位谱lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2); %1024点幅度谱，对数幅度谱，相位谱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: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];%显示连续图形的采样点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点幅度谱，对数幅度谱，相位谱lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2);%1024点幅度谱，对数幅度谱，相位谱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];%显示连续图形的采样点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点幅度谱，对数幅度谱，相位谱lgmagyn=20*log10(abs(ynfft));argyn=angle(ynfft);magyn2=abs(ynfft2);%1024点幅度谱，对数幅度谱，相位谱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];%显示连续图形的采样点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('正态分布噪声'),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