中国地质大学现代数字信号处理报告

现代数字信号处理上机报告 任课老师 姓 名 学 号 专 业 第1题（1）代码：clear all; close all; sample=1000; %取1000个样本点 %误差判定阈值 e1=0.0008; e2=0.002; %N=300;%滤波器阶数 n=1:sample; %产生x,y方向的期望信号 dnx=cos(0.002*pi*n); dny=sin(0.002*pi*n); %产生x,y方向的噪声 xn=sqrt(0.03)*randn(1,sample); yn=sqrt(0.02)*randn(1,sample); %产生x,y方向的信号 x=dnx+xn; y=dny+yn; %产生x方向上观测信号的自相关函数 rxx=xcorr(x); %产生x方向上观测信号与期望信号的互相关函数 rxd=xcorr(x,dnx); %产生y方向上观测信号的自相关函数 ryy=xcorr(y); %产生y方向上观测信号与期望信号的互相关函数 ryd=xcorr(y,dny); %x方向滤波 for Nx=1:10:999 for i=1:Nx for j=1:Nx Rxx(i,j)=rxx(sample-i+j); %x方向信号的自相关矩阵 end end for i=1:Nx Rxd(i)=rxd(sample-1+i); %x方向信号和期望信号的互相关矩阵 end hoptx=inv(Rxx)*Rxd;%由维纳-霍夫方程得到的x方向上的滤波器最优解 ax=1 zeros(1,Nx-1); fx=filter(hoptx,ax,x);%滤波后x方向上的输出 ex=dnx-fx;%误差信号 Ex=sum(abs(ex).2)/sample; %最小均方误差 if(Ex=e1) break; end end %y方向滤波 for Ny=1:10:999 for i=1:Ny for j=1:Ny Ryy(i,j)=ryy(sample-i+j); %y方向信号的自相关矩阵 end end for i=1:Ny Ryd(i)=ryd(sample-1+i); %y方向信号和期望信号的互相关矩阵 end hopty=inv(Ryy)*Ryd;%由维纳-霍夫方程得到的y方向上的滤波器最优解 ay=1 zeros(1,Ny-1); fy=filter(hopty,ay,y);%滤波后y方向上的输出 ey=dny-fy;%误差信号 Ey=sum(abs(ey).2)/sample; %最小均方误差 if(Ey=e2) break; end end%绘制图形%subplot(2,5,1) plot(dnx); axis(0 1000 -1.5 1.5);title(x方向期望信号); subplot(2,5,2) plot(xn); axis(0 1000 -1.5 1.5);title(x方向噪声信号); subplot(2,5,3) plot(x); axis(0 1000 -1.5 1.5);title(x方向观测信号); subplot(2,5,4) plot(fx); axis(0 1000 -1.5 1.5);title(x方向滤波后信号);subplot(2,5,5) plot(ex); axis(0 1000 -1 1);title(x方向误差信号);subplot(2,5,6) plot(dny); axis(0 1000 -1.5 1.5);title(y方向期望信号); subplot(2,5,7) plot(yn); axis(0 1000 -1.5 1.5);title(y方向噪声信号); subplot(2,5,8) plot(y); title(y方向观测信号); subplot(2,5,9) plot(fy); axis(0 1000 -1.5 1.5);title(y方向滤波后信号);subplot(2,5,10) plot(ey); axis(0 1000 -1 1);title(y方向误差信号); figure; plot(dnx,dny,k); hold on; plot(x,y,b:); hold on; plot(fx,fy,g-); axis(-2 2 -2 2)title(最终结果);（2）运行结果：（3）遇到的问题及解决方法：题目要求设计FIR维纳滤波器滤除信号中的高斯噪声，首先需要确定滤波器的阶数。一开始，我直接设置滤波器为300阶。后来经过老师提醒，改成了用误差判定阈值来确定滤波器的阶数。因为滤波器的效果主要是由均方误差来评估的。设计维纳滤波器的过程就是求解维-霍夫方程的过程。观察方程结构，问题就转化成求观测信号的自相关矩阵Rxx、观测信号与期望信号的互相关矩阵Rxd。Matlab自带的xcorr函数可以帮助我们求信号的相关函数。这里需要注意的是，xcorr求得的相关函数是一个偶函数，所以自相关矩阵和互相关矩阵都从相关函数的二分之一处开始取值。求得滤波器的最优解后，我遇到一个问题就是用卷积conv滤波，序列的长度发生了改变，滤波后的输出序列应该截取哪一段？我搜查了一些资料，和同学讨论一下，得出的结论是：截取前面的数据比较好。也就是说，如果输入序列长度为1000，滤波器是200阶，二者卷积之后得到的是1199长度的序列，截取前面的1000个数据就行。因为我们只是对信号的一段数据进行考察，这段数据之后的数据未知，这就导致输出数据的末尾无法精确计算，也就是会产生“边界效应”。从最终结果来看，滤波后点的运动轨迹最后有点偏差，就是因为“边界效应”。filter函数滤波的输出数据和输入数据的长度一样，所以最后我选用filter滤波。第2题（1）代码：clear all;close all;samples=1000;%采样点数N=200;%滤波器阶数Tn=200*pi/3;%信号周期n=0:(Tn*10)/(samples-1):Tn*10;dn=6*sin(0.03*n);%期望信号noise=sqrt(6)*randn(1,samples);%均值为0，方差为6的高斯噪声信号xn=dn+noise;%受干扰的观测信号%绘制期望信号和受干扰观测信号波形图figure(1)plot(dn);hold on;plot(xn,g);title(受干扰信号和期望信号);rxx=xcorr(xn); for i=1:N for j=1:N Rxx(i,j)=rxx(samples-i+j); end end %计算收敛因子u X,B=eig(Rxx); Eig=sort(diag(B); Eigmax=Eig(length(Eig);%观测信号自相关矩阵的最大特征值 u=rand(1)/Eigmax; rdx=xcorr(xn,dn); for i=1:N Rdx(i)=rdx(samples-1+i); endw=zeros(1,N);%加权系数初始化Wn=zeros(samples-N,N);yn=zeros(1,samples);%观测信号初始化yn(1:N)=xn(1:N);en=zeros(1,samples);%误差信号初始化E=zeros(1,samples-N);%迭代800次for i=(N+1):samples x=xn(i-N):i-1); yn(i)=w*x; en(i)=dn(i)-yn(i);%误差信号 E(i)=en(i).2; w=w+2*u*en(i)*x;%Widrow-Hoff LMS算法求权系数 Wn(i-N),:)=w; endfigure(2)plot(dn);hold on;plot(yn,g);title(期望信号和滤波后信号);figure(3)plot(E);title(均方误差变化);figure(4)plot(Wn);title(权系数变化);（2）运行结果：（3）遇到的问题及解决方法：自适应滤波器就是利用前一时刻已经获得的滤波器参数等结果，自动地调整现时刻的滤波器参数以适应信号与噪声未知的或随时间变化的统计特性，从而实现最优滤波。这里我用的是最陡下降法搜索最佳权系数。根据最陡下降法的递推公式，我们要确定收敛因子u和梯度负方向-j。一开始，我直接对收敛因子进行赋值，没有考虑LMS算法的收敛性质。由于u值取得很大，所以得到的结果收敛很快，但跟踪性能很差。后来我根据公式来确定u值，先求得观测信号自相关矩阵的最大特征值max，然后控制u在如下公式范围内。由Widrow-Hoff LMS算法，采用梯度的估计值代替梯度的精确值，最终最陡下降法的递推公式修改为：根据这个公式迭代800次，不断调整权系数。从运行结果可以看出：随着迭代次数的增加，滤波效果越来越好，均方误差逐渐收敛，加权矢量逐渐在横轴波动。第3题（1）代码：clear all;close all;N=1024;%数据长度为Nwn=sqrt(0.36)*randn(1,N);%白噪声wnxn=zeros(1,N);%平稳随机信号xnfor i=1:N-1 xn(i+1)=0.5*xn(i)+wn(i+1);endfigure(1)plot(wn);%绘制白噪声信号图axis(0 N -5 5);xlabel(n);ylabel(wn);title(白噪声);figure(2)plot(xn);%绘制平稳随机信号图axis(0 N -5 5);xlabel(n);ylabel(xn);title(平稳随机信号);%理想功率谱H,W1=freqz(xn);Fs=1024;%采样率f=W1*Fs/(2*pi);figure(3)plot(f,10*log10(abs(H).2);%绘制理想功率谱xlim(0 512)xlabel(f);ylabel(P);title(理想功率谱);x=xn(1:4:N);index=0:round(N/2-1);P1=(abs(fft(x,1024).2/256;Out1=10*log10(P1(index+1);f1=(0:length(Out1)-1)*Fs/length(Out1);figure(4)plot(f1,Out1);xlim(0 512)xlabel(f);ylabel(P);title(256点周期图);%绘制256点周期图P2=(abs(fft(xn,1024).2/1024;Out2=10*log10(P2(index+1);f2=(0:length(Out2)-1)*Fs/length(Out2);figure(5)plot(f2,Out2);xlim(0 512)xlabel(f);ylabel(P);title(1024点周期图);%绘1024点周期图（2）运行结果：（3）遇到的问题及解决方法：用周期图法估计信号的功率谱主要步骤就是对信号求N点的FFT，然后取模的平方再除以N。周期图属于渐近无偏估计，N越大，周期图的统计平均值趋于它的真值。从运行结果来看，1024点的周期图谱估计比256点的周期图谱估计更加接近理想功率谱。第4题（1）代码：clear all;close all;N=1000;%数据长度为1024xn=sqrt(0.36)*randn(1,N);%白噪声%均值为零方差为0.36的白噪声yn=zeros(1,N);zn=zeros(1,N);for i=1:N-1 yn(i+1)=(xn(i+1)+xn(i)/2;endfor j=1:N-1 zn(j+1)=(yn(j+1)-yn(j)/2;endfigure(1)plot(xn);title(xn序列);figure(2)plot(yn);title(yn序列);figure(3)plot(zn);title(zn序列);%zn的均值和方差Ez=mean(zn)Varz=var(zn)Pz,w=pwelch(zn,N);%zn的自功率谱figure(4)plot(w,10*log10(abs(Pz);xlabel(/);ylabel(P/dB);title(zn的自功率谱);（2）运行结果：Ez = -1.2371e-004Varz = 0.0432（3）遇到的问题及解决方法：求zn的均值可以用MATLAB中的mean函数，求方差用var函数。题目给出pwelch函数绘制zn的自功率谱，刚开始我对这个函数不太了解。后来看了MATLAB的函数介绍和网上的一些资料，明白了pwelch的用法。Pwelch的调用形式：Px,w=pwelch(X,WINDOW)，X是信号，WINDOW如果是一个int类型的值，信号X将被分成长度为WINDOW的小段，默认使用长度为WINDOW的汉明窗。最终得到的Px是求得的功率谱，w是归一化频率。第5题（1）代码：clear all;close all;Fs=1;%1HZ采样N=1000;%采样点数index=0:round(N/2-1);whitenoise=sqrt(0.03)*randn(1,1000);%均值为0方差为0.03的高斯白噪声%低通滤波器%f1=0.4;f2=0.47;wp=2*f1/Fs;%截止频率ws=2*f2/Fs;rp=0.8;rs=40;%波纹系数n,wn=cheb1ord(wp,ws,rp,rs);%求切比雪夫1型滤波器最小阶数bz1,az1=cheby1(n,rp,wp);h,w=freqz(bz1,az1,256,Fs);%求解滤波器的频率响应hl=10*log10(abs(h);xn=filter(bz1,az1,whitenoise);%低通滤波a,b=xcorr(xn,unbiased);a=abs(fft(a,N);%高通滤波器%f3=0.1;f4=0.01;wph=2*f3/Fs;%截止频率wsh=2*f4/Fs;rph=0.8;rsh=40;%波纹系数nh,wnh=cheb1ord(wph,wsh,rph,rsh);%求切比雪夫1型滤波器最小阶数bz2,az2=cheby1(n,rph,wph,high);hh,wh=freqz(bz2,az2,256,Fs);%求解滤波器的频率响应hh=10*log10(abs(hh);yn=filter(bz2,az2,whitenoise);%高通滤波c,d=xcorr(yn,unbiased);c=abs(fft(c,N);%绘制图形%fl=linspace(0,Fs,length(xn);fh=linspace(0,Fs,length(yn);figure(1)subplot(2,2,1);plot(w,hl);ylim(-50 5)xlabel(w/Hz);ylabel(h1/dB);title(低通滤波器频谱特性);subplot(2,2,2);plot(fl,10*log10(a);axis(0 0.5 -50 15);xlabel(f/Hz);ylabel(P/dB);title(低通滤波后的功率谱);subplot(2,2,3);plot(wh,hh);ylim(-50 5)xlabel(w/Hz);ylabel(hh/dB);title(高通滤波器频谱特性);subplot(2,2,4);plot(fh,10*log10(c);axis(0 0.5 -50 15);xlabel(f/Hz);ylabel(P/dB);title(高通滤波后的功率谱);figure(2)Pxy,ww=cpsd(xn,yn,N);plot(ww/(2*pi),10*log10(abs(Pxy);xlabel(f/Hz);ylabel(Pxy/dB);title(互功率谱);（2）运行结果：（3）遇到的问题及解决方法：题目要求设计通带在0.1Hz，0.4Hz这一小段范围有重叠的低通滤波器和高通滤波器。我设计的是切比雪夫1型滤波器。在MATLAB中可以用cheb1ord函数求得滤波器的最小阶数，然后用cheby1函数获取滤波器的参数。设计过程中需要解决的问题是如何确定两个滤波器的截止频率。由于这两个滤波器的通带在0.1Hz，0.4Hz重叠，所以可以知道低通滤波器的通带截止频率为0.4Hz，高通滤波器的通带截止平率为0.1Hz。滤波器的阻带截止频率只要遵循奈奎斯特定理即可。滤波器的频谱特性可以用freqz函数求得它的频率响应，然后取幅度谱即可。由于信号的自相关函数和信号的功率谱是一对傅里叶变换，所以求滤波后信号的自功率谱可以先求它的自相关函数再对它进行快速傅里叶变换。从运行结果来看，低通滤波后得到的信号功率谱在0.4Hz以后为0 ，高通滤波后得到的信号功率谱在0.1之前为0。求互功率谱用Pxy,w= cpsd(X，Y，WINDOW)，信号X、Y将被分成长度为WINDOW的小段，默认使用长度为WINDOW的汉明窗。最终得到的Pxy是求得的互功率谱，w是归一化频率。第6题（1）代码：clear all;close all;xn=load(data.mat);N1=6;N2=10;%自相关法Pxx11,w11=pyulear(xn.x,N1);figure(1)subplot(1,2,1);plot(w11,10*log10(abs(Pxx11);xlabel(/);ylabel(P/dB);title(N=6 自相关法估计功率谱);Pxx12,w12=pyulear(xn.x,N2);subplot(1,2,2);plot(w12,10*log10(abs(Pxx12);xlabel(/);ylabel(P/dB);title(N=10 自相关法估计功率谱);%Burg法Pxx21,w21=pburg(xn.x,N1);figure(2)subplot(1,2,1);plot(w21,10*log10(abs(Pxx21);xlabel(/);ylabel(P/dB);title(N=6 Burg法估计功率谱);Pxx22,w22=pburg(xn.x,N2);subplot(1,2,2);plot(w22,10*log10(abs(Pxx22);xlabel(/);ylabel(P/dB);title(N=10 Burg法估计功率谱);%协方差法Pxx31,w31=pcov(xn.x,N1);figure(3)subplot(1,2,1);plot(w31,10*log10(abs(Pxx31);xlabel(/);ylabel(P/dB);title(N=6 协方差法估计功率谱);Pxx32,w32=pcov(xn.x,N2);subplot(1,2,2);plot(w32,10*log10(abs(Pxx32);xlabel(/);ylabel(P/dB);title(N=10 协方差法估计功率谱);%改进的协方差法Pxx41,w41=pmcov(xn.x,N1);figure(4)subplot(1,2,1);plot(w41,10*log10(abs(Pxx41);xlabel(/);ylabel(P/dB);title(N=6 改进的协方差法估计功率谱);Pxx42,w42=pmcov(xn.x,N2);subplot(1,2,2);plot(w42,10*log10(abs(Pxx42);xlabel(/);ylabel(P/dB);title(N=10 改进的协方差法估计功率谱);（2）运行结果：（3）遇到的问题及解决方法：题目要求估计给定数据的AR模型参数。AR模型参数估计方法有4种：自相关法、Burg法、协方差法和修正协方差法。一开始我用的是自相关法，模型阶数k从1开始递推，调用aryule函数求k阶AR模型的参数得到各个系数和白噪声方差Vk。随着模型阶数的增加，求得的白噪声的方差在逐渐收敛。设定一个方差判定阈值V，如果k阶AR模型的白噪声方差小于V就结束递推，输出此刻的模型阶数、系数和白噪声方差。后来老师提醒我，自相关法还是跟经典谱估计一样，默认观测数据以外的值为0，造成分辨率降低。所以我对比了四种AR模型参数估计方法，用pyulear、pburg、pcov和pmcov函数估计数据的功率谱。从运行结果可以看出：自相关法分辨率较低，其他三种方法得到的功率谱分辨率较高；此外，AR模型阶次越高，分辨率越高，但阶次太高，会使估计误差加大造成谱峰分裂。小结通过这次上机实习，我发现自己有很多地方没有做好。在本科阶段，我的数字信号处理基础没学好，而且目前的研究方向跟这个没有太大的联系，所以在做上机题的时候遇到很多问题。有些题虽然根据公式或者PPT上的讲解能够做出来，但是一些深入的问题我还不太了解，还有一些MATLAB的函数用法也不太清楚。经过这次实习，我对现代数字信号处理有了从理论到实际的认识。尤其是维纳滤波器的设计、自适应滤波器的设计、功率谱估计。对MATLAB处理信号也有了进一步了解。这次实习在今后的学习中将会给我带来很多帮助。对课程的建议1.因为很多课程都是最近结课，所以上机时间有点紧。建议每上完一部分内容，在本次课的下一周安排上机，上机报告可以最后一起交。2.课程报告的题目可以早点确定，方便大家早点自学相关内容。