数字信号实验报告.doc

测试信号处理实验报告班 级： 姓 名： 专 业： 学 号： 实验目录1. 有限IIR滤波器的设计11.1 IIR数字滤波器设计原理11.2 步骤及内容11.3 实验结果如下运行结果：21.4 IIR数字滤波器源程序22维纳滤波器32.1 维纳滤波器设计原理32.2. 步骤及内容32.3 实验结果42.4 维纳滤波器源程序63 卡尔曼滤波93.1 卡尔曼滤波设计原理93.2 步骤及内容103.3 实验结果：103.4 卡尔曼滤波程序114. 自适应信号滤波134.1自适应信号滤波器设计原理134.2设计步骤与方法144.3实验结果154.4 自适应滤波器源程序17211. 有限IIR滤波器的设计1.1 IIR数字滤波器设计原理利用双线性变换设计IIR滤波器，首先要设计出满足指标要求的模拟滤波器的传递函数，然后由通过双线性变换可得所要设计的IIR滤波器的系统函数。如果给定的指标为数字滤波器的指标，则首先要转换成模拟滤波器的技术指标，这里主要是边界频率的转换，对指标不作变化。边界频率的转换关系为。接着，按照模拟低通滤波器的技术指标根据相应设计公式求出滤波器的阶数和截止频率；根据阶数查巴特沃斯归一化低通滤波器参数表，得到归一化传输函数；最后，将代入去归一，得到实际的模拟滤波器传输函数。之后，通过双线性变换法转换公式，得到所要设计的IIR滤波器的系统函数。1.2 步骤及内容1) 用双线性变换法设计一个巴特沃斯IIR低通数字滤波器。设计指标参数为：在通带内频率低于时，最大衰减小于；在阻带内频率区间上，最小衰减大于。2) 以为采样间隔，绘制出数字滤波器在频率区间上的幅频响应特性曲线。1.3 实验结果如下运行结果：1.4 IIR数字滤波器源程序clearrp=1;rs=15;wp=.2*pi;ws=.3*pi;wap=tan(wp/2);was=tan(ws/2);n,wn=buttord(wap,was,rp,rs,s);z,p,k=buttap(n);bp,ap=zp2tf(z,p,k);bs,as=lp2lp(bp,ap,wap);bz,az=bilinear(bs,as,.5);h,f=freqz(bz,az,256,1);plot(f,abs(h),r*)title(双线性z变换法获得数字低通滤波器，归一化频率轴);xlabel(omega/2pi);ylabel(低通滤波器的幅频响应);grid;figure;h,f=freqz(bz,az,256,100);ff=2*pi*f/100;absh=abs(h);plot(ff(1:128),absh(1:128),r*);title(双线性z变换法获得数字低通滤波器,频率轴取0,pi/2);xlabel(omega);ylabel(低通滤波器的幅频响应);grid on;2维纳滤波器2.1 维纳滤波器设计原理维纳滤波是一种从噪声中提取信号的最佳线性方法，假定一个随机信号下具有以下形式： （1）其中 表示有用信号，表示噪声。当其输入一个单位脉冲响应为的线性系统时，其输出为： （2）2.2. 步骤及内容实验中，其中， 是零均值方差为的均匀分布白噪声，是与互不相关的均匀分布白噪声，零均值方差为1。维纳最佳滤波器为 （3）冲激响应：则其中为的最佳估值利用 （4）其中L为滤波数据长度。2.3 实验结果1.维纳预测输入样本点个数 4000输入FIR滤波器阶数 10 滤波前的均方误差：0.997589515425905 FIR滤波器的均方误差：0.250057137233902维纳滤波器的均方误差：0.115678278671965通过程序运行结果的均方误差比较，可以看出信号x（n）经过维纳过滤后，噪音已基本消除，与趋于一致。L固定，则N越大，滤波效果越好；N 固定，则L越大，滤波效果越好。2.4 维纳滤波器源程序clear% Program Filter in Matlab Author: huanboma Time: 6/10/2009L = 4000; %Input the number of signal samplesN = 10; %Input FIR Filter ordera = 0.95;K = 40;sigma_a2 = 1-a2;a_ = 1,-a;% Generate random signal x(n)=s(n)+v(n) and the requirement is rou_xx0.03,rou_xs0.01while(1) %s(n) is Source signal, v(n) is Noise, x(n) is the received signalw = sqrt(sigma_a2)*(randn(L,1); %w is the white noise, its variance is sigma_a2v = randn(L,1); %v is the white noise, its variance is 1s(1) = w(1); %x(n) is got through the formula: s(n)=a*s(n-1)+wn(n)for i=1:L-1s(i+1)=a*s(i)+w(i+1);x(i) = s(i)+v(i);endx(L) = s(L)+v(L);%Correlation function r_xx and r_xsfor i=1:K+1rxx(i)=sum(x(i:L).*x(1:L-i+1)/(L-i+1);endr_xx_g=rxx(K+1:-1:2),rxx(1:K+1);for i=1:K+1rxs(i)=sum(x(i:L).*s(1:L-i+1)/(L-i+1);endr_xs_g=rxs(K+1:-1:2),rxs(1:K+1);%Whether or not to meet rou_xx0.03,rou_xs0.01r_xx_t = a.abs(-K:K);r_xx_t(K+1) = r_xx_t(K+1)+1;r_xs_t = a.abs(-K:K);rou_xx=(sum(r_xx_g-r_xx_t).2)/sum(r_xx_t.2);rou_xs=(sum(r_xs_g-r_xs_t).2)/sum(r_xs_t.2);if rou_xx0.03 & rou_xs0.01break;endend % Figure The Theoretical value and Received value of x(n) correlation functionfigure(1);stem(r_xx_t,r);hold onstem(r_xx_g,o,b);legend(theoretical value,received value);title(Theoretical value and Received value of x(n) correlation function);%From the following figure, the theoretical value and received value of x(n)s correlation function are almost the same.%The original signal can achieved through Wiener filter.% Figure The last 100 s(n) and x(n)figure(2);stem(s(L-100:L),r);hold onstem(x(L-100:L),o,b);legend(s(n),x(n);title(The last 100 s(n) and x(n);%From the following figure, there are some differences between the source signal s(n) and received signal x(n).%But they are roughly the same. Because x(n) has noise in same time. % Figure The Theoretical value and Estimated value of h(n)n = 0:N-1;h_t = 0.238*(0.724).n; %h(n)s theoretical valuefor i=1:Nxx(i)=sum(x(i:L).*x(1:L-i+1)/(L-i+1);endfor i=1:NRxx(i,1:N)=xx(i:-1:1),xx(2:N+1-i);endfor i=1:Nxs(i)=sum(x(i:L).*s(1:L-i+1)/(L-i+1);endinvRxx=inv(Rxx);h_g=invRxx*xs; %h(n)s estimated valuefigure(3);stem(h_t,r);hold onstem(h_g,o,b);legend(Theoretical value,Estimated value);title(Theoretical value and Estimated value of h(n);%From the following figure, the value of h(n) can be estimated from Wiener Filter.%And it is similar to the theoretical value of h(n). % Figure The last 100 s(n) and Wiener Filter s_w(n)y_l(1) = s(1); %s_w(n) from Wiener fliterfor i = 1:L-1y_l(i+1) = 0.724*s(i)+0.238*x(i+1);endfigure(4);stem(s(L-100:L),r);hold onstem(y_l(L-100:L),o,b);legend(s(n),Wiener Filter S_l(n);title(The last 100 s(n)(red) and Wiener Filter s_w(n);%From the following figure, the source signal s(n) can be estimated from Wiener Filter, it is called s_w(n).%And it is similar to the original source signal s(n).% Figure The last 100 s(n) and FIR fliter s_f(n)for i=1:9 %s_f(n) from FIR flitery_F(i)=sum(h_g(1:i).*x(i:-1:1);endfor i=10:Ly_F(i)=sum(h_g.*x(i:-1:i-9);endfigure(5);stem(s(L-100:L),r);hold onstem(y_F(L-100:L),o,b);legend(s(n),FIR fliter S_R(n);title(The last 100 s(n) and FIR fliter s_f(n);%From the following figure, the source signal s(n) can be estimated from FIR Filter, it is called s_f(n).%And it is also similar to the original source signal s(n).% Mean Square Errore_x=sum(x-s(1:L).2)/L; %MSE before Filtere_w=sum(y_l-s(1:L).2)/L ; %MSE from Wiener Filtere_f=sum(y_F-s(1:L).2)/L ; %MSE from FIR Filter3 卡尔曼滤波3.1 卡尔曼滤波设计原理卡尔曼滤波是一种递推估计方法，它也是基于最小均方误差准则，但它不需要全部过去的观测数据，只是根据前一个的估计值和最近一个观测数据来估计信号的当前值。在稳态情况下，卡尔曼滤波结果和维纳滤波结果相同。假设已知动态系统的状态方程和量测方程：首先我们要利用系统的过程模型，来预测下一状态的系统。假设现在的系统状态是k，根据系统的模型，可以基于系统的上一状态而预测出现在状态： (1)式(1)中，X(k|k-1)是利用上一状态预测的结果，X(k-1|k-1)是上一状态最优的结果，U(k)为现在状态的控制量，如果没有控制量，它可以为0。到现在为止，结果已经更新，可是，对应于X(k|k-1)还没更新。我们用P表示： (2)式(2)中，P(k|k-1)是X(k|k-1)对应P(k-1|k-1)是X(k-1|k-1)对应的A表示A的转置矩阵，Q是系统过程的式子1，2就是卡尔曼滤波器5个公式当中的前两个，也就是对系统的预测。现在我们有了现在状态的预测结果，然后我们再收集现在状态的测量值。结合预测值和测量值，我们可以得到现在状态(k)的最优化估算值X(k|k)： （3）其中Kg为卡尔曼增益(Kalman Gain)： （4）目前为止，我们已经得到了k状态下最优的估算值X(k|k)。但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束，我们还要更新k状态下X(k|k) （5）其中I 为1的矩阵，对于单模型单测量，I=1。当系统进入k+1状态时，P(k|k)就是式子(2)的P(k-1|k-1)。这样，算法就可以自回归的运算下去。卡尔曼滤波器的原理基本描述了，式子（1-5）5 个基本公式。根据这5个公式，可以很容易的实现计算机的程序。3.2 步骤及内容1. 编制卡尔曼滤波器的通用程序2. 运行卡尔曼滤波器程序，其中U（k）=0，选择L=100，观察并记录实验结果，分析滤波效果3.3 实验结果：3.4 卡尔曼滤波程序clearL=200;Rw=1;Rv=1;xe(1)=0;Pe(1)=1e-12;P(1)=100;for i=1:L A(i)=0.95; B(i)=0; W(i)=1; C(i)=1; V(i)=1; u(i)=0;end r1=randn(1,L);r2=randn(1,L);w=sqrt(Rw)*r1;v=sqrt(Rv)*r2;r=randn(1);x(1)=xe(1)+sqrt(Pe(1);for k=2:L x(k)=0.95*x(k-1)+w(k-1);end y=x+v;for k=1:L yI(k)=C(k)*x(k);end I=eye(1,1);t=1:L; figure(1)plot(t,yI,t,y,:)xlabel(时间 t)title(量测值 yI. y); for k=2:L Q(k-1)=W(k-1)*W(k-1)*Rw; R(k)=V(k)*V(k)*Rv; P1(k)=A(k)*P(k-1)*A(k)+Q(k-1); H(k)=P1(k)*C(k)*inv(C(k)*P1(k)*C(k)+R(k); P(k)=(I-H(k)*C(k)*P1(k); xe(k)=A(k)*xe(k-1)+H(k)*(y(k)-C(k)*A(k)*xe(k-1)+B(k)*u(k-1); ye(k)=C(k)*xe(k);end figure(2)plot(t,P);xlabel(时间 t)title(P)figure(3)plot(t,H);xlabel(时间 t)title(H)figure(4)plot(t,x,t,xe,:);xlabel(时间 t)title(纯信号 x .滤波结果 xe);figure(5)plot(t,yI,t,ye,:)xlabel(时间 t)title(量测值 yI .滤波结果ye);4. 自适应信号滤波4.1 自适应信号滤波器设计原理脉冲响应不变法是实现模拟滤波器数字化的一种直观而常用的方法。他特别适合用于那些对滤波器的时域特性有一定要求的场合。具体的说他可以保证所设计的滤波器的脉冲响应和相应模拟滤波器的冲激响应应在采样点上完全一致。脉冲响应不变法也由此得名。一个模拟滤波器的传递函数可以用有理式表示为 （4-1）通过拉式变换，我们就可以得到他的冲激响应 （4-2）脉冲响应不变法就是要保持脉冲响应不变，即 （4-3）对上式中的冲激响应序列h(n)作z变换，就可以得到数字滤波器的传递函数 （4-4）上式（4-3）和式（4-4），我们已在一些讨论采样序列Z变换的章节中反复接触到。H（Z）就是h（t）理想采样的Z变换。因此，脉冲响应不变的法也成为标准Z变换法。一般来说，h(s)的分母多项式阶次总是大于分子多项式的阶次。假定h(s)没有多重极点（有重极点时要对以下讨论作适当的修正），则式（4-1）可以展成部分分式的形式： （4-5）其中，均为复数， 是 的极点。其拉式变换为 （4-6）由此得到数字滤波器的单位脉冲响应序列 （4-7）数字滤波器的传递函数则为 （4-8）上式经过合并简化，成为一般形式的有理分式传递函数 （4-9）4.2 设计步骤与方法自适应FIR维纳滤波器的LMS算法公式：因此给定初始值h0,每得到一个样本yn,可以得到一组新的滤波器权系数hn。如果信号yn为一个M阶的AR模型，则利用LMS算法可以对AR模型参数进行自适应估计。综上所述，我们可以将脉冲响应不变法设计IIR滤波器的过程归纳为如下几步骤：a) 根据要求，确定数字滤波器和相应模拟滤波器的个临界频率，和衰减，其中 ，T为采样周期。b) 安给定的和设计模拟滤波器传递函数c) 将展成部分分式形式，见式（4-5），求d) 将和代入式（4-8），并化简得到式（4-9），从而完成数字滤波器传递函数H（z）的设计。e) 分析所设计滤波器的频域特性，检查它是否满足要求。上述涉及过程的计算比较复杂，对于一个高阶的滤波器的设计（5阶以上），一般只能借助计算机来完成。本实验要求实验者编制脉冲响应不变数字滤波器设计程序，在计算机上完成上述涉及过程。4.3 实验结果1自适应滤波输入样本点个数 200，输入w的方差 10，输入hI =1.6，输入步长u =0.06输入初始值h0 = 0 可见，u越大，h(n)和s(n)的速度越快。可见a1(n)、a2(n)随着n的增加趋近于a1、a2.4.4 自适应滤波器源程序clearL=input(输入样本点个数 );a=input(输入w的方差 );hI=input(输入hI );u=input(输入步长u );h0=input(输入初始值h0 ); n=1:L;x=randn(1,L);w=sqrt(a)*x;s=hI*x;y=s+w;figure(1)plot(n,s,n,y,:);title(曲线s 曲线y.);ylabel(幅度);xlabel(迭代次数n);hh(1)=0;h(1)=0;%hI+(1-2*u)*(h0-hI);ss(1)=0;sss(1)=0;for i=2:L hh(i)=hh(i-1)+2*u*(s(i-1)-hh(i-1)*x(i-1)*x(i-1); ss(i)=hh(i)*x(i); h(i)=hI+(1-2*u)i)*(h0-hI); sss(i)=h(i)*x(i);end figure(2)plot(n,h,n,hh,:);title(曲线h 曲线hh.);ylabel(幅度);xlabel(迭代次数n);figure(3)plot(n,s,n,ss,:);title(曲线s 曲线ss.);ylabel(幅度);xlabel(迭代次数n);% lms Program start clear all%channel system ordersysorder = 5 ;% Number of system pointsN=2000;inp = randn(N,1);n = randn(N,1);b,a = butter(2,0.25);Gz = tf(b,a,-1);%This function is submitted to make inverse Z-transform (Matlab central file exchange)%The first sysorder weight value%h=ldiv(b,a,sysorder);% if you use ldiv this will give h :filter weights to beh= 0.0976;0.2873;0.3360;0.2210;0.0964;y = lsim(Gz,inp); % Add some noisen = n * std(y)/(10*std(n);d = y + n;totallength=size(d,1);%Take 60 points for trainingN=60 ; % Begin of algorithmw = zeros ( sysorder , 1 ) ;for n = sysorder : N u = inp(n:-1:n-sysorder+1) ;y(n)= w * u;e(n) = d(n) - y(n) ;% Start with big mu for speeding the convergence then slow down to reach the correct weightsif n 20mu=0.32;elsemu=0.15;endw = w + mu * u * e(n) ;end % Check of resultsfor n = N+1 : totallengthu = inp(n:-1:n-sysorder+1) ;y(n) = w * u ;e(n) = d(n) - y(n) ;end hold onplot(d)plot(y,r);title(System output) ;xlabel(Samples)ylabel(True and estimated output)figuresemilogy(abs(e) ;title(Error curve) ;xlabel(Samples)ylabel(Error value)figureplot(h, k+)hold onplot(w, r*)legend(Actual we