现代信号处理例题及matlab代码实现.doc

现代信号处理期末考核作业1 MATLAB仿真均值为0，方差为1的白噪声信号，信号长度N=1024，并用周期图法分别求500、1000和1500次实现的平均功率谱密度，画图。程序代码如下：clear;clear all;N=1024;%数据长度Nfft=1024;%FFT所采用的数据长度n=0:N-1;wn=randn(1,N);%产生随机白噪声subplot(2,2,1);%绘出白噪声序列plot(n,wn);title(白噪声);%500次实现的平均功率谱密度s=zeros(1,N);for i=1:500wn=randn(1,N);%产生随机白噪声Pxx=10*log10(abs(fft(wn,Nfft).2)/N);%Fourier振幅谱平方的平均值，并转换为dbs=s+Pxx;ends=s/500;f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列subplot(222);plot(f,s);xlabel(频率/Hz);ylabel(功率谱/dB);title(500次实现的平均功率谱密度);grid on;%1000次实现的平均功率谱密度s=zeros(1,N);for i=1:1000wn=randn(1,N);%产生随机白噪声Pxx=10*log10(abs(fft(wn,Nfft).2)/N);%Fourier振幅谱平方的平均值，并转换为dbs=s+Pxx;ends=s/1000;f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列subplot(223);plot(f,s);xlabel(频率/Hz);ylabel(功率谱/dB);title(1000次实现的平均功率谱密度);grid on;%500次实现的平均功率谱密度s=zeros(1,N);for i=1:1500wn=randn(1,N);%产生随机白噪声Pxx=10*log10(abs(fft(wn,Nfft).2)/N);%Fourier振幅谱平方的平均值，并转换为dbs=s+Pxx;ends=s/1500;f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列subplot(224);plot(f,s);xlabel(频率/Hz);ylabel(功率谱/dB);title(1500次实现的平均功率谱密度);grid on;实验结果图如下：2仿真如下随机过程：其中：Vn是均值为0，方差为1的Gaussian白噪声过程，为随机相位，在0，2间服从均匀分布。试对其中的正弦波频率进行估计（在不同的数据长度下，N=16，64，128，1024，可使用经典谱估计中的任何一种方法），并讨论数据长度对估计分辨率和平滑特性的影响。解答：使用周期图法对不同数据长度的信号进行估计。程序代码如下：clear all;%*N=16;第一种情况数据长度为16*N = 16;Nfft = 16;n = 0:N-1;xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N);figure(1); subplot(2,1,1);plot(n,xn); ylabel(幅值(V);xlabel(时间(s);title(原始信号);Pxx=10*log10(abs(fft(xn,Nfft).2)/N); Fourier振幅谱平方的平均值，并转换为dbf = 0:length(Pxx)-1; 绘出频率序列subplot(212);plot(f,Pxx);xlabel(频率/Hz);ylabel(功率谱/dB);title(周期图N=16);grid on;%*N=64*N = 64;Nfft = 64;n = 0:N-1;xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N);figure(2); subplot(2,1,1);plot(n,xn); ylabel(幅值(V);xlabel(时间(s);title(原始信号);Pxx=10*log10(abs(fft(xn,Nfft).2)/N);f = 0:length(Pxx)-1;subplot(212);plot(f,Pxx);xlabel(频率/Hz);ylabel(功率谱/dB);title(周期图N=64);grid on;%*N=64*N = 128;Nfft = 128;n = 0:N-1;xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N);figure(3); subplot(2,1,1);plot(n,xn); ylabel(幅值(V);xlabel(时间(s);title(原始信号);Pxx = 10*log10(abs(fft(xn,Nfft).2)/N);f = 0:length(Pxx)-1;subplot(212);plot(f,Pxx);xlabel(频率/Hz);ylabel(功率谱/dB);title(周期图N=128);grid on;%*N=64*N = 1024;Nfft =1024;n = 0:N-1;xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N);figure(4); subplot(2,1,1);plot(n,xn); ylabel(幅值(V);xlabel(时间(s);title(原始信号);Pxx = 10*log10(abs(fft(xn,Nfft).2)/N);f = 0:length(Pxx)-1;subplot(212);plot(f,Pxx);xlabel(频率/Hz);ylabel(功率谱/dB);title(周期图N=1024);grid on;实验结果如下：数据长度N=16时的谱估计结果数据长度N=64时的谱估计结果数据长度N=128时的谱估计结果数据长度N=1024时的谱估计结果3计算并画图描绘下列函数离散时间傅立叶变换的幅度（幅频）特性：n=0，131n=其它 n=0，131n=其它 频率范围2，2解答：程序代码如下：clear all;close all;N=32;dt=1;%设置最大点数n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;%离散Fourier变换矩阵t=n*dt;%*xn信号xn=exp(j*n*pi/5);figure(1);subplot(311);plot(t,xn);title(xn信号);Xk=xn*WNnk;%对xn进行Fourier变换magXk=abs(Xk);phaXk=angle(Xk);k=0:length(magXk)-1;subplot(312);plot(k/(N*dt),magXk*2/N);title(xn信号的振幅谱);subplot(313);plot(k/(N*dt),unwrap(phaXk);ylabel(相位角/rad);title(xn信号的相位谱);%*yn信号yn=sin(n*pi/5);figure(2);subplot(311);plot(t,yn);title(yn信号);Yk=yn*WNnk;%对xn进行Fourier变换magYk=abs(Yk);phaYk=angle(Yk);k=0:length(magYk)-1;subplot(312);plot(k/(N*dt),magYk*2/N);title(yn信号的振幅谱);subplot(313);plot(k/(N*dt),unwrap(phaYk);ylabel(相位角/rad);title(yn信号的相位谱);实验结果图如下：信号的离散时间傅里叶变换信号的离散时间傅里叶变换4设Rxx(0)=7.24，设Rxx(1)=3.6，试确定如下一阶MA的参数值用两种方法： 直接用MA方程： 谱分解方法。最后描绘用MA模型得到的谱估计。解答：MA方程方法：由，c0=7.24, c1=3.6得b0=1.3379，b1=2.6907。程序代码如下：b=1.3379 2.6907; %MA系统系数w=linspace(0,pi,512);H=freqz(b, w);%产生信号的频域响应Ps=abs(H).2;%计算得到功率谱plot(w./(2*pi),Ps);xlabel(频率/Hz);ylabel(功率谱/dB);title(直接用MA方程方法得到的谱估计);谱分解方法：程序代码如下：clear all;N=456;B1=1 2; A1=1;w=linspace(0,pi,512);%采用自协方差法对AR模型参数进行估计%y1=filter(B1,A1,randn(1,N).*zeros(1,200),ones(1,256);Py11,F=pcov(y1,1,512,1); %AR(1)的估计%Py13,F=periodogram(y1,512,1);%*MA模型*%y=zeros(1,256);for i=1:256y(i)=y1(200+i);endny=0:255;z=fliplr(y);nz=-fliplr(ny); nb=ny(1)+nz(1);ne=ny(length(y)+nz(length(z);n=nb:ne;Ry=conv(y,z);R1=zeros(2,2);r1=zeros(2,1);for i=1:2r1(i,1)=-Ry(260+i);for j=1:2R1(i,j)=Ry(260+i-j);endendR1; r1; a1=inv(R1*R1)*R1*r1; %利用最小二乘法得到的估计参数%对MA的参数进行估计%A1;A11=1,a1; %AR的参数的估计值B11=fliplr(conv(fliplr(B1),fliplr(A11); %MA模型的分子y21=filter(B11,A1,randn(1,N);%.*zeros(1,200),ones(1,256)