页数：125页文件大小：2.82 MB来源： MATLAB实验教程ppt课件
本课件配合数字信号处理课程，用12个实例演示脉冲响应不变法与双线性变换法设计IIR带通、带阻滤波器，FIR滤波器设计、模拟滤波器原型与巴特沃斯滤波器，双音多频检测，以及自适应滤波LMS、维纳滤波与卡尔曼滤波。
MATLAB中的信号处理工具箱提供了很多数字信号处理中需要用到的函数及解决方法。在数字信号处理中所学习到的问题，诸如滤波器的设计、自适应滤波，维纳滤波、卡尔曼滤波等理论都可以通过MATLAB仿真得到实现和验证。本章提供了12个MATLAB在数字信号处理课程学习中的应用实例，可作为学习数字信号处理课程的参考。
知识架构
1．实验目的
T为采样周期。若以Ha(s)及H(z)分别表示ha(t)的拉氏变换及h(n)的Z变换，即：根据采样序列Z变换与模拟信号拉氏变换的关系得:上式表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S平面到Z平面的变换，正是以前讨论的拉氏变换到Z变换的标准变换关系。
脉冲响应不变法特别适用于用部分分式表达的传递函数，模拟滤波器的传递函数若只有单阶极点，且分母的阶数高于分子阶数，则可表达为部分分式形式：其拉氏反变换为：对ha(t)采样就得到数字滤波器的单位脉冲响应序列：再对h(n)取Z变换，得到数字滤波器的传递函数：
3．仿真思路在MATLAB中，可以用下列函数辅助设计IIR数字滤波器。
clearwp=6*pi*10^3;ws=9*pi*10^3;ap=1,as=15;Fs=30*10^3;wp1=wp/Fs;ws1=ws/Fs;[N,WC]=cheb1ord(wp,ws,ap,as,'s');[b,a]=cheby1(N,ap,WC,‘s’)；[bz,az]=impinvar(b,a,Fs)；w0=[wp1,ws1]；Hx=freqz(bz,az,w0)；[H,W]=freqz(bz,az)；dbHx=-20*log10(abs(Hx)/max(abs(H)))；plot(W,abs(H));xlabel('相对频率');ylabel('幅频');grid1．实验目的
查巴特沃斯归一化低通滤波参数表可得：设计巴特沃斯低通滤波器。阶数N计算如下。数字滤波器系统函数H(z)为：
A=0.09036B1=1.2686, C1=−0.7051B2=1.0106, C2=−0.3583B3=0.9044, C3=−0.2155可见H(z)是由3个二阶滤波器H1(z)、H2(z)、H3(z)级联组成的，如图4.3所示。图4.3滤波器H(z)的组成图2)带阻IIR数字滤波器实现过程4．程序代码1)心电图信号滤波代码x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,...-4,-4,-6,-6,-2, 6,12, 8, 0,-16,...-38,-60,-84,-90,-66,-32,-4,-2,-4, 8,...12,12,10,6,6,6,4,0,0,0,...0,0,-2,-4,0,0,0,-2,-2,0,...0,-2,-2,-2,-2,0];%存有高频干扰心电图信号序列k=1;close all;Figure(1)subplot(2,2,1)axis([0 56-100 50]);n=0:55;stem(n,x,'.');hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('心电图信号采样序列x(n)');B=[0.09036 2*0.09036 0.09036];A=[1.2686-0.7051];A1=[1.0106-0.3583];A2=[0.9044-0.2155];while(k<=3)y=filter(B,A,x);%对信号x进行滤波x=y;if k==2A=A1;endif k==3A=A2;endk=k+1;endsubplot(2,2,3)axis([0 56-15 5]);n=0:55;stem(n,y,'.');hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('y(n)');title('三级滤波后的心电图信号');2)带阻滤波器设计程序代码wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000;wp=[wp1,wp2]/(Fs/2);%利用奈奎斯特频率归一化ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');%求滤波器阶数[num,den]=ellip(N,rp,rs,wc,'stop');%求滤波器传递函数[H,W]=freqz(num,den);%频率响应plot(W*Fs/(2*pi),abs(H));grid;%绘制网格xlabel('频率/Hz');ylabel('幅值')5．运行结果与分析从图4.4中可以看出低通滤波器除高频噪声的滤波效果。由该图可以看出，滤波前的心电图信号很不整齐，滤波后的心电图信号变化很有规律，消除了毛刺等噪声的影响。图4.5中的带阻滤波器从通带到阻带变化很陡峭，滤波效果不错。6．思考题
1．实验目的
假设理想频响的采样是H(k)，(k=0,1,…,N-1)则其IDFT为：则FIR滤波器的系统函数可写为：所以当采样点数N已知后，只要采样值H(k)确定，则系统函数H(z)就可以确定，要求的FIR滤波器就设计出来了。2)线性相位FIR滤波器的约束条件若分别是对幅度函数H(ω)
和相位函数φ(ψ)的第k个抽样点。在ω∈[0，2π]内的N个样点，则约束条件为:该程序除了给出设计的滤波器的理想低通频谱取样点图、频率取样法设计的低通脉冲响应图以及3．仿真思路低通FIR滤波器的MATLAB设计技术指标为：一线性相位低通滤波器中截止频率为15Hz，τ=12s。MATLAB中滤波器设计的函数采用的是归一化频率，首先要将截止频率化为归一化频率：
滤波器的低通衰耗特性图外，还给出一个由两个正弦频率混合得到的信号通过该滤波器后得到的时域波形和幅频波形图。同时为了增大阻带衰耗，采用增加过渡带样点，如在两边过渡带取一个值为0.4的过渡带取点。4．程序代码%线性相位FIR滤波器的频率取样法的设计Fs=50;%采样频率N=25;%取样点数tao =(N-1)/2;wc=2*15/50 ;%N是奇数h(n)为偶数对称,采用N为奇数的Ⅰ型,取样化为归一化频率
m=fix(wc/(2/N)+1) ;%取整omega=[0:N-1]*2*pi/N ;%取N个频率样点ABS_H=[ones(1,m),0.4,zeros(1,N-2*m-1),0.4,ones(1,m-1)];H=ABS_H.*exp (-j*tao.*omega);h=ifft(H);figure(1);%第1幅图subplot(2,1,1);stem(omega,abs(H));%画出低通频谱采样点图title('理想低通频谱采样点')subplot(2,1,2);stem([0:N-1],real(h));%画出低通脉冲响应图title('频率采样法设计的低通脉冲响应')Hh=fft(h,512);%对信号进行傅里叶变换Hh_db=20*log10(abs(Hh));figure(2);%第2幅图omega=[0:255]*pi/256;plot(omega,Hh_db(1:256));title('滤波器的低通衰耗特性')set(gca,'ytick',[(-120:20:-20),-3,0]);t=(0:99)/Fs;s=sin(2*pi*t*10)+sin(2*pi*t*20);%给出频率分别是10Hz和20Hz的混合信号进行测试sf=fft(s,512);%对信号进行傅里叶变换y=conv(s,h);%相乘运算yf=fft(y,512);%对信号进行傅里叶变换figure(3);%第3幅图subplot(3,1,1);plot(0.5.*([1:length(y)]-1),y);%画出输入的混合信号时域波形图title('输入的混合信号波形图')subplot(3,1,2);w1=(0:255)/256*(50/2);plot(w1,abs(sf(1:256)));%画出输入的混合信号的幅频图title('混合信号的幅频波形图')subplot(3,1,3);plot(w1,abs(yf(1:256)));title('混合信号通过滤波器后的幅频波形图')5．运行结果与分析理想低通频谱采样点如图4.6所示。理想低通滤波器的衰耗特性如图4.7所示。由图4.8知，测试频率是10Hz和20Hz，通过滤波器后，20Hz大于截止频率，所以被滤除掉。而10Hz小于截止频率所以可通过滤波器。通过上三图可看出，滤波器的各项时域和频域指标均达到了要求。图4.6理想低通频谱采样点图4.7理想低通滤波器的衰耗特性
图4.8输入的混合信号波形图6．思考题
1．实验目的
3)带通切比雪夫Ｉ型滤波器的设计对于带通滤波器的设计可以通过低通滤波器进行频率变换来实现。下面介绍确定带通滤波器的技术指标的步骤。
4) MATLAB中设计带通切比雪夫Ⅰ型滤波器的频率变换函数
3．仿真思路本实验中给出3个正弦信号，频率分别为300Hz、500Hz、1000Hz，将这3个信号进行叠加，显示合成信号频谱，再通过低通巴特沃斯滤波器分离出300Hz和500Hz的信号，再通过带通切比雪夫Ｉ型滤波器识别出1000Hz的信号，并显示低通和带通滤波器的图形及经过滤波器后的信号频谱图，然后经过扬声器播放合成音、通过低通后和通过带通后的声音。4．程序代码t=1/10000:1/10000:3;%设置3个频率的正弦信号s1=sin(2*pi*300*t);s2=sin(2*pi*500*t);s3=2*sin(2*pi*1000*t);
s=s1+s2+s3;% 3个正弦信号叠加figure(1);subplot(3,1,1); ss=fft(s,4096);SS=(abs(ss(1:1:2049)));%求合成信号频谱k1=0:2048;w1=(1/4096)*k1*10000;%取0...Fs/2的部分plot(w1,SS); grid%画频谱图title('求原信号频谱');% % %通过低通滤波器ws1=1000;wp1=600;wc=5000;wp=wp1/wc;ws=ws1/wc;[n,wn]=buttord(wp,ws,1,30);%巴特沃斯低通滤波[b,a]=butter(n,wn);subplot(3,1,2);sb=3*filter(b,a,s);%合成信号通过低通滤波器ssb=fft(sb,4096);SSb=abs(ssb(1:1:2049));%求频谱k1=0:2048;w1=(1/4096)*k1*10000;%画频谱图plot(w1,SSb);gridtitle('经过低通滤波器后的信号频谱');% % %通过带通滤波器ws1=[600 1400];%设计一个阻带为600~1400Hzwp1=[980 1020];%通带为980~1020Hz的带通滤波器wc=5000;wp=wp1/wc;ws=ws1/wc;[n,wn]=cheb1ord(wp,ws,4,45,'s');%切比雪夫I型滤波器[bd,ad]=cheby1(n,4,wn,'bandpass');subplot(3,1,3);sd=3*filter(bd,ad,s);%合成信号通过带通滤波器ssd=fft(sd,4096);SSd=abs(ssd(1:1:2049));%求频谱k1=0:2048;w1=(1/4096)*k1*10000;plot(w1,SSd);grid%;画频谱图title('经过带通滤波器后的信号频谱');figure(2)%画低通滤波器的幅频、相频图freqz(b,a,4096,10000);axis([0,4000,-70,1]);figure(3)%画带通滤波器的幅频、相频图freqz(bd,ad,4096,10000);axis([0,4000,-70,1]);sound(s,10000);pause(5);sound(sb,10000);pause(5);sound(sd,10000);5．运行结果与分析低通滤波器幅频响应及相频响应如图4.9所示。带通滤波器幅频响应及相频响应如图4.10所示。经过滤波器的信号频谱如图4.11所示。图4.9图4.10图4.11经过滤波器的信号频谱
通过以上对滤波器的设计，可以得出通过设计所需信号的滤波器能滤去其余干扰信号，而且滤波器还可以用来对信号进行变换、调制、解调、去噪、识别信号等。从程序代码中可以看出，只要写入正确的参数，就可以利用MATLAB中的函数去设计所需的滤波器，进而运用于实际中。6．思考题低通滤波器分离合成信号的原理是什么？
1．实验目的
这里的特性接近矩形，是指通带频率响应段与过渡频率响应段的夹角接近直角，通常为钝角，如果该角为直角，则为理想滤波器。所谓的滤波器的零点就是将该点的值代入传递函数后，传递函数的值为零。所谓的函数的极点就是将该点的函数值代入传递函数后，传递函数的值为无穷大。滤波器的增益是传递函数表达式前的常数。若系统的传递函数表示为：则滤波器的零点为[Z
3．仿真思路首先绘制Butterworth低通模拟原型滤波器的幅频平方响应曲线，阶数分别为2、5、10、20，然后再设计一个Butterworth带通滤波器，技术指标为：W1=2000Hz，W2=3000Hz，两侧的过渡带宽为500Hz，带通波纹为Rp=1dB，阻带衰减为Rs=60dB，阶数为N=14。Butterworth函数调用格式为：[z，p，k]=buttap(n)。其中n为滤波器的阶数，z、p、k分别为滤波器的零点极点和增益。传递函数具有下面形式：
滤波器没有零点，极点为[P
switch icase 1,N=2;case 2;N=5;case 3;N=10;case 4;N=20;end[z,p,k]=buttap(N);[b,a]=zp2tf(z,p,k)[H,w]=freqs(b,a,n);magH2=(abs(H)).^2;hold on;plot(w,magH2);endxlabel('w/wc');ylabel('|H(jw)|^2');title('butterworth模拟原型滤波器');text(1.5,0.18,'n=2');text(1.3,0.08,'n=5');text(1.16,0.08,'n=10');text(0.93,0.98,'n=20');grid on;2)程序2wp=[2000 3000]*2*pi;%通带临界频率ws=[1500 3500]*2*pi;%阻带临界频率Rp=1;%通带最大衰减Rs=60;%阻带最小衰减[N,Wn]=buttord(wp,ws,Rp,Rs,‘s’)；%计算阶数和截止频率[b,a]=butter(N,Wn,'s');w=linspace(1,4000,2000)*2*pi;H=freqs(b,a,w);plot(w/(2*pi),20*log10(abs(H)));xlabel('Frequency(Hz)');ylabel('Magnitude(dB)');title(['Butterworth Analog bandpass filter order N=',num2str(N)]);grid on;5．运行结果与分析由图4.12和图4.13可知，滤波器的幅度平方特性随着频率单调下降；随着滤波器的阶数增大，其幅度特性也渐渐接近于矩形，符合设计要求的滤波器。6．思考题巴特沃斯滤波器具有什么特点？图4.12图4.131．实验目的
模拟滤波器按幅度特性可分为低通滤波器、高通滤波器、带通滤波器和带阻滤波器。2)模拟滤波器设计方法模拟滤波器设计方法已经相当成熟，并且有若干典型的模拟滤波器可供选择，如巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器和贝塞尔滤波器等。这些滤波器的特点各异，巴特沃斯滤波器具有单调下降的幅频特性；切比雪夫滤波器的幅频特性在通带和阻带内有波动，可以提高选择性；贝塞尔滤波器在通带内有较好的线性相位特性；椭圆滤波器的选择性相对于前3种滤波器是最好的。
在MATLAB中，提供了上述所有滤波器的设计函数，从而大大降低了滤波器的设计难度。设计模拟滤波器主要有两种方法。一种是首先调用buttap、cheblap等函数设计模拟低通滤波器原型，然后通过lp2lp、lp2hp等函数进行频率变换，将模拟低通滤波器原型转换为高通、带通和带阻等滤波器。另一种是根据滤波器性能参数，首先调用buttord、cheb2ord等阶数估计函数，对满足技术指标的滤波器的最低阶数进行估计，然后再调用滤波器完全设计函数进行模拟滤波器设计。2)信号传输的条件
所谓信号无失真传输是指输入信号通过系统后，输出信号的幅值和输入信号的幅值成正比。允许有一定的延时，但没有波形上的畸变。因此，系统的频率响应H(jw)满足下面的特性。即信号无失真传输的条件是：系统的幅频响应|H(jw)|应为常数，相频响应∠ H(jw)应与ω成比。3．仿真思路巴特沃斯模拟低通滤波器的平方幅频响应函数为：
在有关模拟滤波器设计的MATLAB程序中，经常遇到一些特定函数。[b，a]=zp2tf(z，p，k)为将模拟原型滤波器函数(如buttap)设计出的零点z、极点p和增益k形式转换为传递函数(transfer function)形式。通过[H，ω ]=freqs(b，a，ω )求出传递函数形式表示的滤波器对应于频率点ω的复数频率响应H(包括实部和虚部)，这里ω为一个向量，表示对应的角频率。若该函数不写输出变量，则执行后绘出该滤波器的幅频响应和相频响应图。4．程序代码
n=0:0.01:2;for i=1:4;switch icase 1;N=3case 2;N=5;case 3N=10;case 4N=15;end;[z,p,k]=buttap(N);[b,a]=zp2tf(z,p,k);[h,w]=freqs(b,a,n);magh=abs(h);subplot(2,2,i);plot(w,magh);axis([0 2 0 1]);xlabel('w/wc');ylabel('|H(jw)|^2');title(['Butterworth analog filter N=',num2str(N)]);grid on;end5．运行结果与分析由图4.15可以看出，巴特沃斯滤波器的频率特性曲线无论在通带还是在阻带内，都是频率的单调函数。因此，当通带边界处满足指标要求，通带内有余量，可以将精度均匀地分布在整个通带内。为此，选择传递函数具有等波纹特性的滤波器，就能在阶数较低的情况下满足系统要求。图4.15巴特沃斯模拟低通滤波器的幅频特性曲线6．思考题ButterWorth滤波器的边缘是否陡峭对滤波器的性能有何影响？仿真中，边缘是否陡峭与什么参数有关？
1．实验目的
它的行频和列频唯一标识。标准频率的选定不是任选的，而是设计时考虑了拨号音、忙音及交流电源的杂音，这些杂音均不在上述DTMF信号的范围内。同时为了避免标准频率之间的谐波干扰出现，故所列最高标准频率小于最低标准频率的3次谐波。DTMF是用两个特定的单音频组合信号来代表数字信号以实现其功能的一种编码技术。两个单音频的频率不同，代表的数字或实现的功能也不同。这种电话机中通常有16个按键，其中有10个数字键0～9和6个功能键*、#、A、B、C、D。由于按照组合原理，一般应有8种不同的单音频信号。因此可
采用的频率也有8种，故称之为多频。根据CCITT的建议，国际上采用的频率为697Hz、770Hz、852Hz、941Hz、1209Hz、1336Hz、1477Hz和1633Hz等8种。用这8种频率可形成16种不同的组合，从而代表16种不同的数字或功能键，具体组合如图4.16所示。图4.16 DTMF编码
解码部分用滤波器组实现。使用低通及高通滤波器来完成。每个数字都包含低频和高频两个频率分量，分别设计低通和高通滤波器来分离。通过检测出的有效频率，查找表，可以得到解码数字，输出该数字。3．仿真思路若滤波器的输入、输出都是离散时间信号，那么该滤波器的单位冲激响应h(n)也必然是离散的，这种滤波器称为数字滤波器。当用硬件实现时，所需的元件是乘法器、延时器和相加器；而用MATLAB软件实现时，仅仅需要线性卷积程序即可。
数字滤波器的作用是对输入信号进行滤波，就如同信号通过系统一样。此处使用低通及高通滤波器来完成。每个数字都包含低频和高频两个频率分量，分别设计低通和高通滤波器来分离。低频率分量697Hz、770Hz、852Hz和941Hz的最小间隔为73Hz，设计低通滤波器，使滤波器的通带截止频率介于每两个频率中间，使得检测某一频率时其他频率分量可以很好地被隔离。在这里低通滤波器设计4个，在滤波时逐个调用从低到高一旦检测到有效分量立即终止，保证不出现误检。高频率分量为1209Hz、1336Hz、1477Hz和1633Hz，最小间
隔为127Hz，类似于低通，设计4个高通滤波器，检测时则是从高到低检测，一旦检测到有效分量即终止。分别滤波后，根据循环次数来确定对应的数字。4．程序代码%DTMF信号代表的16个数tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68];f1=[697,770,852,941];%行频率向量f2=[1209,1336,1477,1633];%列频率向量TN=input('键入3位电话号码=') ; %输入3位数字TNr=0;%接收端电话号码初值为0m=3;%分别对每位号码数字处理：产生信号，发声，检测for m=1:3;d=fix(TN/10^(3-m)); %计算出第m位号码数字TN=TN-d*10^(3-m);forp=1:4;forq=1:4;if tm(p,q)==abs(d); break, end%检测与m位号码相符的列号qendif tm(p,q)==abs(d); break, end%检测与m位号码相符的列号pendn=0:400; %为了发声，加长序列x=sin(2*pi*n*f1(p)/8000)+sin(2*pi*n*f2(q)/8000);subplot(3,3,m);plot(n,x);grid;title('未滤波的波形');sound(x,8000); %发出声音pause(0.2) %相邻号码响声之间加0.2秒停顿for k=1:4;wpz=[0.175 0.2 0.225 0.2375];%低频的4个低通截止数字频率wsz=0.3;rp=3;rs=60;[N,wc]=buttord(wpz(k),wsz,rp,rs);%低通数字滤波器[Bz,Az]=butter(N,wc,'low');y=filter(Bz,Az,x);%滤波后的单频信号z=max(abs(y));%滤波后信号在离散域所取的最大值subplot(3,3,m+3),plot(n,y);;grid;title('低通后的波形');xlabel('n');ylabel('y')ifz>0.7,break,end%查找行号endfor l=1:4;wpz2=[0.4 0.3625 0.325 0.3];%低频的4个高通截止数字频率wsz2=0.25;rp2=3;rs2=60;[N2,wc2]=buttord(wpz2(l),wsz2,rp2,rs2);%高通数字滤波器[Bz2,Az2]=butter(N2,wc2,'high');y2=filter(Bz2,Az2,x);%滤波后的单频信号z2=max(abs(y2));%滤波后信号在离散域所取的最大值subplot(3,3,m+6);plot(n,y2);;grid;title('高通后的波形');xlabel('n');ylabel('y')if z2>0.7,break,end;%查找列号endTNr=TNr+tm(k,(5-l))*10^(3-m);%将3位电话号码表示成一个3位数enddisp('接受断检测到的号码为：')disp(TNr)5．运行结果与分析从图4.17中可知，3个数字的低频部分频谱相同，高频部分不同，跟假设匹配。通过对滤出频率分析，查找表，可以得出输出的数字。图4.17数字‘1’、‘2’、‘3’的时域和滤波后的波形6．思考题双音多频滤波器常用来解决什么问题？
1．实验目的
到级联形式的转换以及二者的比较。直接形式用延迟元件和乘法器以及加法器以给定的形式直接实现差分方程，设M=N=4，则差分方程为：在这种形式中将系统函数H(z)写成实系数二阶子系统的乘积形式。首先把分子和分母多项式的根解出，然后把每一对共轭复根或任意两个实根组合在一起，得到二阶子系统。在以下原理中假设N为偶
本程序由1个主函数和3个子函数组成，主函数完成由下面的差分方程描述的直接型滤波器到级联型滤波器的转换。15y(n)+10y(n−1)+ 3y(n−2)−5y(n−3)−2y(n−4)=x(n)−4x(n−1)+10x(n−2)−25x(n−3) +16x(n−4)子函数dir2cas实现将直接形式转化为级联形式，子函数stepseq(n0,n1,n2)，实现阶跃函数u(n-n0)，n1≤n＜n2。子函数casfilter是编制实现级联形式滤波器的函数。级联结构中每一个二阶网络决定一对零点和一对极点，可以通过灵活地调整对应零、极点的系数来
改变一对零、极点的位置，相对于直接型结构，其优点是调整方便，此外，级联结构中后面的网络输出不会流到前面，运算误差的积累相对于直接型也小。4．程序代码1)主程序b=[1-4 10-25 16];a=[15 10 3-5-2];[b0,B,A]=dir2cas(b,a)stepin=stepseq(0,0,5);format long;hcas=casfilter(b0,B,A,stepin)hdir=filter(b,a,stepin)format;figure(1)subplot(2,1,1)stem(hcas); hold onplot([0 6],[0 0]);axis([0 6-1.8 1]);title('级联型的阶跃响应'); hold off;subplot(2,1,2)stem(hdir);hold on;plot([0 6],[0,0]);axis([0 6-1.8 1]);title('直接型的阶跃响应'); hold off;2) stepseq函数function [x,n]=stepseq(n0,n1,n2)%Generate x(n)=u(n-n0);n1<=n<=n2%[x,n]=stepseq(n0,n1,n2)STEP = 0.01;n = n1:STEP:n2;x = n > n0;3) casfilter函数function y=casfilter(b0,B,A,x)[K,L]=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=x;for i=1:Kw(i+1,:)=filter(B(i,:),A(i,:),w(i,:)); endy=b0*w(K+1,:);4) dir2cas函数function [b0,B,A]=dir2cas(b,a)b0=b(1);b=b0/b0;a0=a(1);a=a/a0;b0=b0/a0;M=length(b);N=length(a);if N>Mb=[b zeros(1,N-M)];else if N<Ma=[a zeros(1,N-M)];elseNM=0;endK=floor(N/2);B=zeros(K,3);A=zeros(K,3);ifK*2==Mb=[b 0];a=[a 0];endbroots=cplxpair(roots(b));aroots=cplxpair(roots(a));Fori=1:2:2*K;Brow=broots(i:1:i+1,:);Brow=real(poly(Brow));B(fix((i+1)/2),:)=Brow;Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow));A(fix((i+1)/2),:)=Brow;End5．运行结果与分析b0 =0.0667B =1 0 01 0 0A =1 0 01 0 0由图4.18可知，这两种形式结构的滤波器本质上是完全相同的，但相比于直接型结构，级联形式的滤波器可以更灵活地调整零、极点的系数来改变一对零、极点的位置，调整方便，而且级联结构中后面的网络输出不会流到前面，运算误差的积累也很小，因而得到了广泛应用。
图4.18级联型滤波器和直接型滤波器的阶跃响应6．思考题
1．实验目的
图4.19最小均方误差准则的自适应横向滤波器原理图该自适应滤波器的输入矢量为：加权矢量为为：滤波器的输出为：y(n)相对于滤波器期望输出d(n)的误差为：
根据最小均方误差准则，最佳的滤波器参量应使得性能函数——均方误差为最小。假定输入信号x(n)和期望响应d(n)是联合平稳过程，那么在时刻n的均方误差是加权矢量的二次函数，其表示式为：均方误差是权向量的二次函数，它是一个上凹的抛物面，具有唯一的最小值，调解权向量使得均方误差最小，相当于沿抛物面下降寻找最小值。可以用梯度法来求该最小值，对权向量W求导得到均方误差的梯度为：
在性能曲面上最佳权矢量对应点的梯度等于零，即：该方程称为正则方程，由此解出最佳权向量称为维纳解。利用上式求解，需要精确地知道输入信号和期望信号的先验统计知识，而且还要对矩阵求逆运算。最陡下降法可避免求逆运算，它通过递推的方式寻求加权矢量的最优值，是LMS算法的理论基础。首先设置一个W的初值W
最小均方算法是一种很有用很简单的估计梯度的方法，其突出特点是计算量小、易于实现，且不要求脱线计算。LMS最核心的算法是使用平方误差代替均方误差，即所以将上式代入最陡下降法迭代计算权矢量的公式得：则LMS算法的基本关系式为：
3．仿真思路设计一个二阶加权自适应横向滤波器，对一个正弦信号进行滤波。实验通过设置不同的收敛因子μ，由MATLAB程序图形观察滤波效果。讨论μ的重要性。4．程序代码clear all ;fs =10000;t=0：1/fs：1;sn =sin(2 * pi * t);%产生初始信号n= randn (size (t));%产生高斯噪声xn = sn + n;%信号加噪声w = [0 0.5];%设置权初值u=0.00026;%设置收敛因子for i =1:length(t)-1for i =1:length(t)-1y (i+1)=n(i:i+1)*w';%噪声通过滤波器输出为ye (i+1)=xn (i+1)-y(i+1);w = w+2*u*e(i+1)*n(i:i+1);%权的变化公式endsubplot (3,1,1)plot (t, xn )%输出信号加噪声图形title ('带噪声原始信号')grid；subplot (3,1,2)plot (t , sn)title('原始正弦信号')subplot(3,1,3)plot (t , e)title('滤波结果')5．运行结果与分析收敛因子μ决定收敛速度及稳定性。图4.20设置的收敛因子为0.00026，图4.21设置的收敛因子为0.0026，图4.22设置的收敛因子为0.026。由图可见，收敛因子的选择对滤波器的性能有很大影响。图4.20，滤波效果不错，但是收敛较慢；图4.21，收敛速度加快，但滤波效果不如图4.20；图4.22收敛速度很快，但滤波效果太差。可见，选择合适的收敛因子对于滤波性能有很大影响。收敛因子的选择，其实就是在滤波性能及速度之间折中。6．思考题
图4.20 μ =0.00026图4.21 μ =0.0026图4.22 μ=0.026时正弦加噪信号的滤波
1．实验目的
为：由上式可知，频率域采样点数N必须大于等于时域离散信号的长度M(即N≥M)，才能使时域不产生混叠，则N点的离散傅里叶反变换得到的序列就是原序列项x(n)。对比上面叙述的频域采样定理，得到一个有用的结论：频域采样，时域信号周期延拓。2)频域采样定理的验证给定信号如下：
编写程序分别对频谱函数在区间[0，2π]上等间隔采样32点和16点，得到：3．仿真思路直接调用MATLAB函数fft计算再分别对进行32点和16点IFFT，得到：分别画出的幅度谱，并绘图显示的波形，进行对比和分析。
就得到频谱函数在[0，2π]上的32点和16点频率采样。4．程序代码M=26;N=32;n=0:26;n1=0:M/2;n2=M/2+1:M;xa=n1+1;xb=27-n2;x=[xa,xb];subplot(3,2,2);stem(n,x);title('三角波序列x(n) ')w=linspace(0,pi,1000);%设定频率向量X=x*exp(-j*n'*w);subplot(3,2,1);plot(w/pi,abs(X));title('序列x(n)的连续幅度谱')X32k=fft(x,32);k1=[0:31];subplot(3,2,3);stem(k1,abs(X32k));axis([0,15,0,200]);title('32点频率采样幅度谱')x32n=ifft(X32k);nx32=[0:31];subplot(3,2,4);stem(nx32,x32n);title('32点IDFT')X16k=fft(x,16);k2=[0:15];subplot(3,2,5);stem(k2,abs(X16k));axis([0,8,0,200]);title('16点频率采样幅度谱')x16n=ifft(X16k);nx16=[0:15];subplot(3,2,6);stem(nx16,x16n);axis([0,20,0,28]);title('16点IDFT')5．运行结果与分析由图4.23可知，在一定的条件下，可以由频域离图4.23对不同信号的采样
1．实验目的
可能精确地重现或对信号做出尽可能精确的估计，而对所伴随噪声进行最大限度的抑制。维纳滤波器就是这种滤波器的典型代表之一。1)维纳滤波概述维纳(Wiener)是用来解决从噪声中提取信号的一种滤波方法。这种线性滤波问题可以看成是一种估计问题或一种线性估计问题。如果一个线性系统的单位样本响应为h(n)，当输入一个随机信号x(n)，且其中x(n)表示信号，v(n)表示噪声，则输出y(n)为：
希望x(n)通过线性系统h(n)后得到的y(n)尽量接近于s(n)，因此称y(n)为s(n)的估计值，即：则维纳滤波器的输入-输出关系可用图4.24表示。图4.24维纳滤波器的输入-输出关系用h(n)进行过滤问题实际上是一种统计估计问题。一般地，从当前的和过去的观察值x(n)，x(n −1)，x(n − 2)，…估计当前的信号值称为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值称为外推或预测；从过去的观察值，估计过去的
信号值称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。用e(n)表示信号的真实值和估计值之间的误差，显然e(n)可能是正值，也可能是负值，并且它是一个随机变量。因此用它的均方误差来表达误差是合理的，均方误差最小即它的平方的统计期望最小。采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。2)维纳-霍夫方程的求解
为了按最小均方误差准则来确定维纳滤波器的冲激响应h(n)，令ξ(n)对h(j)的导数等于零，即可得维纳-霍夫方程：其中：维纳-霍夫方程右端的求和范围没有具体标明，实际上有3种情况。
③因果无限冲激响应(因果IIR)维纳滤波器，i从0到+∞取正整数值。上述3种情况下维纳-霍夫方程的解法不同，这里只描述FIR维纳滤波器的求解。设滤波器冲激响应序列的长度为N，冲激响应矢量为：滤波器输入数据矢量为：则滤波器的输出为：这样，上述所示的维纳滤波器的标准方程可写成：
其中，P是s(n)与x(n)的互相关函数，它是一个N维列矢量；R是x(n)的自相关函数，是N阶方阵。利用求逆矩阵的方法直接求解上式，得：opt表示“最佳”，这就是FIR维纳滤波器的冲激响应。3．仿真思路假设一个点目标在x、y平面上绕单位圆做圆周运动，由于外界干扰，其运动轨迹发生了偏移。其中，x方向的干扰为均值为0、方差为0.05的高斯噪声；y方向的干扰为均值为0、方差为0.06的高斯噪声。
clear;clf;sita=0:pi/249.5:2*pi;xnoise=sqrt(0.05)*randn(1,500);%产生x轴方向噪声ynoise=sqrt(0.06)*randn(1,500);%产生y轴方向噪声x=cos(sita)+xnoise;%产生x轴方向观测信号y=sin(sita)+ynoise;%产生y轴方向观测信号%产生维纳滤波中x方向上观测信号的自相关矩阵rxx=xcorr(x);for i=1:100for j=1:100mrxx(i,j)=rxx(500-i+j);endendclear;clf;sita=0:pi/249.5:2*pi;xnoise=sqrt(0.05)*randn(1,500);%产生x轴方向噪声ynoise=sqrt(0.06)*randn(1,500);%产生y轴方向噪声x=cos(sita)+xnoise;%产生x轴方向观测信号y=sin(sita)+ynoise;%产生y轴方向观测信号%产生维纳滤波中x方向上观测信号的自相关矩阵rxx=xcorr(x);for i=1:100for j=1:100mrxx(i,j)=rxx(500-i+j);endendxd=cos(sita);%产生维纳滤波中x方向上观测信号与期望信号的互相关矩阵rxd=xcorr(x,xd);for i=1:100mrxd(i)=rxd(499+i);endhoptx=inv(mrxx)*mrxd‘;%由维纳-霍夫方程得到的x、方向上的滤波器最优解fx=conv(x,hoptx);%滤波后x方向上的输出nx=sum(abs(xd).^2);eminx=nx-mrxd*hoptx;%x方向上最小均方误差%产生维纳滤波中y方向上观测信号的自相关矩阵ryy=xcorr(y);for i=1:100for j=1:100mryy(i,j)=ryy(500-i+j);endendyd=sin(sita);%产生维纳滤波中y方向上观测信号与期望信号的互相关矩阵ryd=xcorr(y,yd);for i=1:100mryd(i)=ryd(499+i);endhopty=inv(mryy)*mryd‘;%由维纳-霍夫方程得到的y方向上的滤波器最优解fy=conv(y,hopty);%滤波后y方向上的输出ny=sum(abs(yd).^2);eminy=ny-mryd*hopty;%y方向上最小均方误差subplot(2,4,1)plot(xd);title('x方向期望信号');subplot(2,4,2)plot(xnoise);title('x方向噪声信号');subplot(2,4,3)plot(x);title('x方向观测信号');subplot(2,4,4)n=0:500;plot(n,eminx);title('x方向最小均方误差');subplot(2,4,5)plot(yd);title('y方向期望信号');subplot(2,4,6)plot(ynoise);title('y方向噪声信号');subplot(2,4,7)plot(y);title('y方向观测信号');subplot(2,4,8)plot(n,eminy);title('y方向最小均方误差');figure;plot(xd,yd,'k');hold on;plot(x,y,'b:');hold on;plot(fx,fy,'g-');title('最终结果');5．运行结果与分析x轴和y轴信号、噪声及最小均方误差如图4.25所示。滤波后得到的信号与原始信号和噪声信号的对比如图4.26所示，滤波后的结果与期望信号还是很接近的，整体上达到了最优滤波的效果。6．思考题
图4.25信号、噪声及最小均方误差滤
图4.26最终结果
1．实验目的
卡尔曼滤波原理首先要引入一个离散控制过程的系统，该系统可用一个线性随机微分方程来描述：X(k)=AX(k−1)+BU(k)+W(k)再加上系统的测量值：Z(k)=H X(k)+V(k)上两式子中，X(k)是k时刻的系统状态，U(k)是k时刻对系统的控制量。A和B是系统参数，对于多模型系统，它们为矩阵。Z(k)是k时刻的测量值，H是测量系统的参数，对于多测量系统，H为矩阵。W(k)和V(k)分别表示过程和测量的噪声。
它们被假设成高斯白噪声，其协方差分别是Q、R这里假设它们不随系统状态变化而变化。由于满足上面的条件(线性随机微分系统，过程和测量都是高斯白噪声)，卡尔曼滤波器是最优的信息处理器。下面来估算系统的最优化输出。首先利用系统的过程模型预测下一个状态的系统。假设现在的系统状态是k，根据系统的模型，可以基于系统的上一状态而预测出现在状态：X(k|k-1)=AX(k-1|k-1)+BU(k)上式中，X(k|k-1)是利用上一个状态预测的结果，X(k-1|k-1)是上一个状态最优的结果，U(k)为现在状
态的控制量，如果没有控制量，它可以为0。到现在为止，系统结果已经更新了，可是对应于X(k|k-1)的协方差还没更新。用P表示协方差：P(k|k−1)=A P(k−1|k−1) A'+Q上式中，P(k|k−1)是X(k|k−1)对应的协方差，P(k−1|k−1)是X(k−1|k−1)对应的协方差，A‘表示A的转置矩阵，Q是系统过程的协方差。上面两个公式就是卡尔曼滤波器5个公式当中的前两个，也就是对系统的预测。有了现在状态的预测结果，再收集现在状态的测量值。结合预测值和测量值，可以得到现在状态(k)
的最优化估算值X(k|k)：X(k|k)=X(k|k−1)+Kg(k)(Z(k)−H X(k|k−1))其中Kg为卡尔曼增益(Kalman Gain)：Kg(k)= P(k|k−1) H'/(H P(k|k−1) H' + R)至此为止，已经得到了k状态下最优的估算值X(k|k)。但是，为了要令卡尔曼滤波器不断地运行下去，直到系统过程结束，还要更新k状态下X(k|k)的covariance：P(k|k)=(I-Kg(k)H)P(k|k-1)其中I为1的矩阵，对于单模型单测量，I=1。根据上述5个公式，可以很容易地实现计算机程序。3．仿真思路把房间看成是一个系统，然后对这个系统建模。已知房间的温度跟前一时刻的温度相同，所以A=1。没有控制量，所以U(k)=0。因此得出：X(k|k-1)=X(k-1|k-1)P(k|k-1)=P(k-1|k-1) +Q因为测量的值是温度计的，跟温度直接对应，所以H=1。X(k|k)= X(k|k-1)+Kg(k) (Z(k)-X(k|k-1))Kg(k)= P(k|k-1) / (P(k|k-1) + R)P(k|k)=(1-Kg(k))P(k|k-1)现在模拟一组测量值作为输入。假设房间的真实温度为25度，模拟200个测量值，这些测量值的平均值为25度，但是加入了标准偏差为几度的高斯白噪声(在图4.27中为蓝线)。为了令卡尔曼滤波器开始工作，需要告诉卡尔曼两个零时刻的初始值X(0|0)和P(0|0)。它们的值不用太在意，随便给一个就可以了，因为随着卡尔曼的工作，X会逐渐收敛。但是对于P，一般不要取0，因为这样可能会令卡尔曼完全相信给定的X(0|0)是系统最优的，从而使算法不能收敛。本例选了X(0|0)=1度，P(0|0)=10。
该系统的真实温度为25度，图4.27中用蓝线表示。图4.27中红线是卡尔曼滤波器输出的最优化结果(该结果在算法中设置了Q=1e−6，R=1e−1)。4．程序代码clearN=200;w(1)=0;w=randn(1,N)x(1)=0;a=1;for k=2:N;x(k)=a*x(k-1)+w(k-1);endV=randn(1,N);q1=std(V);Rvv=q1.^2;q2=std(x);Rxx=q2.^2;q3=std(w);Rww=q3.^2;c=0.2;Y=c*x+V;p(1)=0;s(1)=0;For t=2:N;p1(t)=a.^2*p(t-1)+Rww;b(t)=c*p1(t)/(c.^2*p1(t)+Rvv);图4.27卡尔曼滤波结果
本章小结本章通过12个实例讲述了MATLAB在数字信号处理中的应用。通过本章的学习，将会对MATLAB的数字信号处理工具箱有一定了解，并能作为进一步学习数字信号处理理论的基础。
Thanks书山有路勤为径，学海无涯苦作舟。好好学习，天天向上
下载地址：{{downurl}}
{{{a}}}
{{#b}}{{{b}}}
{{/b}} {{/title}} {{#nolog}}登陆注册
{{/nolog}} {{#djsx}}点我刷新
{{/djsx}} {{#gocz}} {{/gocz}} {{#ffqr}} {{/ffqr}}