Matlab在数字信号处理中的运用.ppt

《Matlab在数字信号处理中的运用.ppt》由会员分享，可在线阅读，更多相关《Matlab在数字信号处理中的运用.ppt（68页珍藏版）》请在三一办公上搜索。

1、第七讲 Matlab在数字信号 处理中的运用,数字信号处理（DSP）：是指利用计算机或专用处理设备，以数值计算的方法对信号进行采集、变换、综合、估值与识别等加工处理，借以达到提取信息和便于应用的目的。Matlab是数字信号处理技术实现的重要手段。通过Matlab，可以进行数字信号处理问题理论上的分析和算法开发；配合Simulink，用户可以进行数字信号处理系统的设计和仿真。,本章内容：,离散时间信号与系统数字滤波器概述IIR滤波器的设计,信号,连续信号 指时间连续、幅度连续的信号，又称模拟信号，数学上表示为一个时间连续函数f(t);离散信号 指时间离散，幅度仍然可以连续的信号。可由f(t)时域

2、取样得到，数学上表示为一个时间离散函数f(n);,数字信号,指时间离散而且幅度也离散的信号，可对模拟信号处理获得：模拟信号 取样、量化、编码 数字信号 上述过程又称为脉冲编码调制，这是一个典型的AD变换过程；计算机处理和存储的信号全部是数字信号，通常需要将数字信号还原为模拟信号，过程为：数字信号 解码、反量化、重建 模拟信号 这个过程又称为 DA变换。,1.离散时间信号与系统,1）离散信号及其Matlab实现单位抽样序列(单位冲激)(n)：仅在n=0时取值为1.在Matlab中，产生N点的单位抽样序列，我们利用zeros函数来实现：x=zeros(1,N);x(1)=1;,例7-1 产生一个3

3、2点的，右移20各单位的冲击序列。程序如下：clear all;N=32;k=20;x=zeros(N);x(k)=1;xn=0:N-1;stem(xn,x);,单位阶跃序列(n),Matlab实现：x=ones(1,N);右移m，则：正弦序列正弦序列定义：x(n)=Asin(2fnTs+)Matlab实现：n=0:N-1;x=A*sin(2*pi*f*n*Ts+fai);,例7-3：产生一个频率为150Hz，幅度为0.45，初始相位为35度的正弦波，信号持续时间为5s。,Fs=2000;t=1/Fs:1/Fs:5;f=150;A=0.45;Fai=35/180*pi;X=A*sin(2*pi

4、*f.*t+Fai);plot(t(1:100),X(1:100);xlabel(time(sec);ylabel(sin2pi ft);title(150Hz sin wave);disp(按任意键开始播出5秒的150Hz正弦波.);pause;,sound(X,Fs);disp(播放结束，下面将音频数据存盘为C:my50HzSIN.wav);wavwrite(X,Fs,C:my50HzSIN.wav);clear;R,Rs=wavread(C:my50HzSIN.wav);sound(R,Rs);,Matlab实现：n=0:N-1;x=a.n;例：x(n)=实现程序：n=0:10;x=(0

5、.9).n;stem(n,x),实指数序列：x(n)=,随机序列,由于许多实际的序列并不能用数学式来描述，我们将这些序列成为随机序列，并用相应的概率密度函数来表示。Matlab提供了两个随机数产生函数：rand(1,N)产生0,1上均匀分布的随机序列，N为长度；randn(1,N)产生均值为0，方差为1的高斯 随机序列，即白噪声序列。,例7-4 产生随机序列,%generate random sequencen=200;xn1=rand(1,n);xn2=randn(1,n);subplot(2,2,1),stem(xn1);xlabel(n);ylabel(x(n);title(rand);

6、gridsubplot(2,2,2);hist(xn1,10);xlabel(n);ylabel(x(n);title(均匀分布的概率密度);gridsubplot(2,2,3),stem(xn2);xlabel(n);ylabel(x(n);title(randn);gridsubplot(2,2,4),hist(xn2,10);xlabel(n);ylabel(x(n);title(高斯分布的概率密度)grid,2）波形发生器,Matlab内部提供了大量的函数用以产生噪声及常用的信号波形，这些信号在信号处理中非常重要。方波函数squaresquare函数有两种调用格式：x=square(t

7、);x=square(t,duty);square函数产生周期为2*pi，幅度为1的方波。duty为占空比，即信号为正值的区域在一个周期内所占的百分比。,例如：x=square(2*pi*f*t+fai,duty)%t为时间取样序列%f为方波的基频率%fai为方波的初始相位%duty为占空比，为0-100之间的数。例如 duty=20，即占空比为20%x 返回的幅度为1的矩形波样值序列,三角函数sawtooth：,x=sawtooth(t);x=sawtooth(t,width);例如：x=sawtooth(2*pi*f*t+fai,width)%f为三角波的基频率%fai为三角波的初相位%w

8、idth为宽度，为0-1之间取值的尺度参数。若width=0.3，表示在三角波的一个周期内上升沿占30%；若width=0.5，产生对称的三角波；若width=1，产生锯齿波。%返回的幅度为1的锯齿波样值序列,其它波形发生函数：,chirp线性调频信号dirichlet周期信号sinc傅立叶反变换rectpuls非周期的、单位高度的矩形信号gauspuls高斯调制正弦脉冲tripuls三角形脉冲信号pulstrain高斯调制正弦脉冲、非周期矩形脉冲 和非周期三角脉冲voc压控振荡信号,3）序列的操作,信号相加：是一种对应的样本与样本之间的相加，表示为：x1(n)+x2(n)=x1(n)+x2(

9、n)用运算符“+”实现，要求参与运算的两个序列长度必须相等。信号相乘：对应采样值之间的相乘，表示为：x1(n)*x2(n)=x1(n)x2(n)用数组运算符“.*”实现，要求参与运算的两个序列长度必须相等。,倍率：每一个采样值乘以一个常数a，表示为ax(n)=ax(n)用算术运算符“*”实现。折叠：x(n)的每个样本都对n=0翻转，得到一个折叠后的序列y(n)。y(n)=x(-n)由flghr(x)实现。,样本和：将n1和n2之间的所有样本x(n)加起来：由sum(x(n1):x(n2)函数实现。样本积：将n1和n2之间的所有样本x(n)乘起来：由prod x(n1:n2)函数实现。,4）线性

10、系统及其Matlab实现,线性系统的基本概念：连续时间线性时不变系统的表示：拉普拉斯变换描述的传递函数形式；一阶微分方程组描述的状态空间形式离散时间线性系统的表示：Z变换描述的传递函数形式；一阶差分方程组描述的状态空间形式,时域响应工具箱函数,当系统由传递函数、状态方程给出时，Matlab 给出了专门求解系统单位冲激相应的函数：impulse(sys):计算连续系统的冲激响应 dimpulse:计算离散系统的冲激响应,例如：计算 的冲激响应。sys=tf(1,1,1,0.5)Transfer function:1-s2+s+0.5impulse(sys);,离散系统模型时域表示,filter函

11、数：利用递归滤波器或非递归滤波器对数据进行滤波。因为一个离散系统可以看作是一个滤波器，系统的输出就是输入经过滤波器滤波的结果。y=filter(b,a,x);表示由向量b和a组成的系统对输入x进行滤波，系统的输出为y；y,zf=filter(b,a,x,zi);zi表示输出信号的初始状态，zf表示该函数返回的系统的最终状态向量。,impz函数：直接给出系统的单位冲激相应，调用格式：impz(b,a),例7-5：当系统的输入差分方程为：y(n)-0.8y(n-1)-0.5y(n-2)=0.7x(n)+0.3(n-1),分别利用filter函数和impz函数求系统的单位冲激相应。clear all

12、;pulse=1,zeros(1,63);b=0.7,0.3;a=1,-0.8,-0.5;h1=filter(b,a,pulse);h2=impz(b,a,64);subplot(2,1,1),stem(h1),title(filter fuction);subplot(2,1,2),stem(h2),title(impz fuction);,传递函数响应,freqs函数：用于计算并画出连续系统的幅频响应和相频响应。常用的调用格式为：h=freqs(b,a)h=freqs(b,a,w)其中：b为传递函数H(s)分子多项式系数，a为分母多项式系统，w是指定计算频率点序列，如果w省略，则自动取20

13、0个频率点作计算。h为返回值，是对应于频率点序列w的复频率响应。,freqz函数：用于计算并画出离散系统的幅频响应和相频响应。该函数使用基于FFT算法计算系统传递函数响应模型中的系数向量a和b。常用调用格式为：h,f=freqz(b,a,n,fs)其中：b为传递函数H(z)分子多项式系数，a为分母多项式系数，n为指定计算频率点数(由于采用FFT算法，n常取 2的幂次方，以提高计算速度)，fs为离散系统的采样频率，h为对应于频率点序列f的复频率响应。如无输出变量，则自动作出幅频响应和相频响应图，f为记录频率点数。,2.数字滤波器概述,1）数字滤波器的数学描述和分类 数字滤波器(DF，Digita

14、l Filter)在数字信号处理中起着重要的作用。在信号的处理，检测与参数的估计方面，数字滤波器是使用最为广泛的一种系统。,由于信号通常夹杂噪声及无用信号成分，所以必须将这些干扰成分滤除。滤波器可以对信号进行筛选，只让特定信号通过。一般而言，噪声信号往往是高频信号。而经典滤波器正是假定有用信号与噪声具有不同的频段，所以利用经典的滤波器可以将噪声滤除。但是如果信号和噪声频谱相互重叠，那么利用经典的滤波器就不能发挥作用。现代滤波器的作用是从含有噪声的数据记录中估计除信号的某些特征或信号本身，那么估计出来的信号和源信号相比，就具有更高的信噪比。,滤波是信号处理的基础，滤波运算是信号处理中的基本运算。

15、因此，滤波器的实际问题也是数字信号处理中的主要问题。数字滤波器的本质是将一组输入的数字序列通过一定的运算后，转变成为另一组输出的数字序列。,滤波器的数学描述有两种方法：,差分方程：系统函数：,滤波器的分类按计算方法：递归系统，非递归系统按冲激响应长度：IIR，FIR按频带分类：低通、高通、带通、带阻,2）设计步骤,按照实际需要确定滤波器的性能要求。用一个因果稳定的系统函数（传递函数）去逼近这个性能要求，这种传递函数可以分为两类：IIR和FIR。用一个有限精度的运算去实现这个传递函数。包括选择运算结构：如级联型、并联型、卷积型等。选择合适字长和有效的数字处理方法。,3.IIR滤波器的设计,一个N

16、阶IIR滤波器的传递函数可以表示为：传递函数的设计就是确定系数ai、bi或零、极点ci、di，以使滤波器满足给定的性能要求。,利用模拟滤波器的理论：先设计一个合适的模拟滤波器，然后变换成满足预定指标的数字滤波器。由于模拟的网络综合理论已经发展的很成熟，已经产生了许多高效率的设计方法。很多常用的模拟滤波器不仅有简单而严格的设计公式，而且其参数已经表格化，设计起来非常方便、准确。因而可以将这些理论继承下来，作为设计数字滤波器的工具。,1）设计方法一般有两种：,最优化设计方法 确定一种最优准则，如最小均方误差准则，使设计出来的实际频率响应的幅频特性的均方误差最小。在此最佳准则下，通过迭代运算求滤波器

17、的系数ai，bi。这里，第一种方法用的比较普遍，而随着计算机技术的发展，最优化设计方法的使用逐渐增多。,IIR滤波器设计一般要借助于模拟滤波器原型，Matlab中设计IIR滤波器的典型步骤如下：先根据数字滤波器的设计要求，获得相应的模拟低通滤波器的技术指标；进行模拟低通滤波器原形设计；在S域上进行频率转换得到相应的带通（高通、带阻等）模拟滤波器；最后再进行一个从S域到Z域的映射变换，就得到了所需的数字滤波器。Signal Processing Toolbox对这几步都提供有实现的 函数。,1）模拟滤波器设计,根据滤波器滤波功能的不同，可以将滤波器分为：低通滤波器高通滤波器带通滤波器带阻滤波器

18、一个离散时间系统可以用它的频率响应 来描述，Signal Processing Toolbox提供有一个函数freqz用于计算数字滤波器的频率响应：,实际中理想滤波器是不可能实现的，在通带和阻带之间应该存在一个过渡带，而且通带和阻带内滤波器的幅频响应也很难做到严格的1或0，应该设定一个门限。实际滤波器设计往往采用 和 分别表示通带截止频率和阻带截止频率，和 分别表示通带允许的最大衰减和阻带应达到的最小衰减，其关系式为：,滤波器的6个重要的通带、阻带参数是:通带截至频率（Hz）;:阻带起始频率（Hz）:通带内波动（dB），即通带内所允许的最大 衰减;:阻带内最小衰减（dB）设采样率为，那么可将以

19、上参数中的频率参 数转换为归一化角频率参数：,:通带截至角频率（rad/s）:阻带起始角频率（rad/s）通过这些参数就可以进行离散滤波器的设计了。,归一化数字频率：如果数字信号是对模拟信号以周期T采样得到，那么数字频率 和相应模拟频率 的关系为：。称为归一化的数字频率。将归一化数字频率乘上 就得到数字频率，单位为：；如果乘以采样速率的一半，得到的是一个Hz表示的频率。,例7-6：某滤波器的传递函数为：利用Matlab画出它的幅频响应和相频响应的曲线。,编写M文件,b=0.0181 0.0543 0.0543 0.0181;a=1.000-1.7600 1.1829-0.2781;h w=fr

20、eqz(b,a);subplot(2,1,1);plot(w/pi,abs(h);xlabel(归一化频率(pi rad/sample);ylabel(幅频响应);subplot(2,1,2);plot(w/pi,angle(h)/pi);xlabel(归一化频率(pi rad/sample);ylabel(相频相应);,幅频相频响应,2）巴特沃斯模拟滤波器设计,模拟滤波器设计方法分为三类：Butterworth z,p,k=buttap(n)chebyshew、型 z,p,k=cheblap(n,Rp)z,p,k=cheb2ap(n,Rp)椭圆滤波器 z,p,k=ellipap(n,Rp,R

21、s),阶数的确定：在模拟滤波器原型设计函数之前，首先需要根据设计要求确定合适的滤波器阶数，阶数取得越高滤波器的性能越好，对于实现来说越不经济。如何选择满足要求的最小阶数，就是一个重要问题。,利用Matlab提供的巴特沃斯滤波器设计函数”buttord”即可求出所需要的滤波器最小阶数和3dB截止频率。命令格式为：n,Wn=buttord(Wp,Ws,Rp,Rs)其中：Wp是通带截止频率 Ws是阻带起始频率 Rp是通带内波动（dB）Rs是阻带内最小衰减(dB)n是巴特沃斯滤波器最低阶数 Wn是巴特沃斯滤波器截止频率,巴特沃斯滤波器系数计算 由巴特沃斯滤波器的阶数n以及截止频率Wn可以计算出对应传递

22、函数H(z)的分子分母系数，Matlab提供的命令是：巴特沃斯低通滤波器系数计算：b,a=butter(n,Wn)巴特沃斯高通滤波器系数计算 b,a=butter(n,Wn,high),巴特沃斯带通滤波器系数计算 b,a=butter(n,W1,W2)巴特沃斯带阻滤波器系数计算 b,a=butter(ceil(n/2),W1,W2,stop)其中：n为低通滤波器阶数 Wn为低通滤波器截止频率 b为 H(z)的分子多项式系数 a为H(z)分母多项式系数,例7-7 巴特沃斯滤波器设计实例,1.采样率为8000Hz,要求设计一个低通滤波器程序如下：f_N=8000;f_p=2100;f_s=2500

23、;R_p=3;R_s=25;Ws=f_s/(f_N/2);Wp=f_p/(f_N/2);n,Wn=buttord(Wp,Ws,R_p,R_s);b,a=butter(n,Wn);freqz(b,a,1000,8000);subplot(2,1,1);axis(0 4000-30 3),2.采样率为8000Hz,要求设计一个高通滤波器，程序如下：f_N=8000;f_p=1000;f_s=700;R_p=3;R_s=20;Ws=f_s/(f_N/2);Wp=f_p/(f_N/2);n,Wn=buttord(Wp,Ws,R_p,R_s);b,a=butter(n,Wn,high);freqz(b,

24、a,1000,8000);subplot(2,1,1);axis(0 4000-30 3),3.采样率为10000Hz,要求设计一个带通滤波器，程序如下：f_N=10000;f_p=1000,1500;f_s=600,1900;R_p=3;R_s=20;Ws=f_s/(f_N/2);Wp=f_p/(f_N/2);n,Wn=buttord(Wp,Ws,R_p,R_s);b,a=butter(n,Wn);freqz(b,a,1000,10000);subplot(2,1,1);axis(0 5000-30 3),4.采样率为10000Hz,要求设计一个带阻滤波器，程序如下：f_N=10000;f_

25、p=1000,1500;f_s=1200,1300;R_p=3;R_s=20;Ws=f_s/(f_N/2);Wp=f_p/(f_N/2);n,Wn=buttord(Wp,Ws,R_p,R_s);b,a=butter(n,Wn,stop);freqz(b,a,1000,10000);subplot(2,1,1);axis(0 5000-35 3),3）频率变换,选定模拟滤波器的类型和确定模拟滤波器的阶数后，就可以进行模拟低通滤波器原型的设计了。但还需要在s平面上进行原型变换，得到相应符合设计要求的模拟滤波器，为下一步转换为数字滤波器作准备。,滤波器原型到低通模拟滤波器：lp2lp 调用格式如下：

26、bt,at=lp2lp(b,a,Wo)其中，参数Wo指定截至频率，系统可以由传递函数来描述（b，a分别代表分子系数向量和分母系数向量）滤波器原型到高通模拟滤波器：lp2hp滤波器原型到带通模拟滤波器：lp2bp滤波器原型到带阻模拟滤波器：lp2bs,例7-8：设计一模拟低通Butterworth滤波器要求：通带截止频率：阻带截止频率：带内最大衰减：带内最小衰减：,clear;Wp=5000*2*pi;Ws=10000*2*pi;Rp=3;Rs=30;n,Wn=buttord(Wp,Ws,Rp,Rs,s);z,p,k=buttap(n);b0,a0=zp2tf(z,p,k);b,a=lp2lp(

27、b0,a0,Wn);h,w=freqs(b,a);plot(w/(2*pi),20*log10(abs(h);grid on;xlabel(频率(Hz);ylabel(幅频响应(dB);xlim(0,12000);ylim(-40,1);,4）从模拟滤波器到数字滤波器,从模拟滤波器变换到数字滤波器，也就是要从模拟滤波器的传递函数H(s)得到数字滤波器的传递函数H(z)，实际是从s平面到z平面的映射变换。常用的方法有2个：,冲激响应不变法 bz,az=impinvar(b,a,fs)bz,az=impinvar(b,a)双线性变换法 zd,pd,kd=bilinear(z,p,k,fs)numd

28、,dend=bininear(num,den,fs)Ad,Bd,Cd,Dd=bilinear(A,B,C,D,fs),例7-9：用冲激响应不变法和双线性变换法将上题设计出的模拟滤波器转换为数字滤波器，采样频率为40kHz.,冲击响应不变法：Wp=5000*2*pi;Ws=10000*2*pi;Fs=40000;Rp=3;Rs=30;n,Wn=buttord(Wp,Ws,Rp,Rs,s);z,p,k=buttap(n);b0,a0=zp2tf(z,p,k);b,a=lp2lp(b0,a0,Wn);bz,az=impinvar(b,a,Fs);freqz(bz,az);,双线性变换法clear;Wp=5000*2*pi;Ws=10000*2*pi;Fs=40000;Rp=3;Rs=30;n,Wn=buttord(Wp,Ws,Rp,Rs,s);z,p,k=buttap(n);b0,a0=zp2tf(z,p,k);b,a=lp2lp(b0,a0,Wn);bz,az=bilinear(b,a,Fs);freqz(bz,az);,