MATLAB数字信号处理ppt课件

第5章 使用MATLAB实现数字信号 处理 本章主要内容如下： 51 数字信号处理基本内容及相应的MATLAB 工具 52 信号通过系统的时域分析 53 信号通过系统的频域和Z域分析 54 滤波器设计 55 频谱分析 5.1 数字信号处理基本内容 及相应的MATLAB工具数字信号处理的基本内容通常分为两部分： 离散时间信号与系统分析主要涉及离散时间信号与系统的时域、频域表示，以及信号通过系统的时域、频域分析及其变换域分析。MATLAB函数库中提供了filter, conv, convmtx, fft ,ifft,freqz, impz, zplane等与之相应的函数。 数字滤波器设计和谱分析 数字滤波器设计包括了无限冲激响应(IIR)和有限冲激 响应（FIR）滤波器设计，谱分析又可进一步分为线性 谱分析和非线性谱分析。MATLAB为此提供了多种成 熟算法的相应函数以及极为丰富的设计工具。5.2 时域分析卷积，滤波，单位冲激响应5.2.1 卷积 MATLAB提供 conv函数实现标准的一维信号卷积 : 例如，若系统h(n)为 >>h=1 1 1 输入序列x(n)为 >>x=1 1 1 则x(n)经过系统h(n)后的MATLAB实现为： >>conv(h,x) 或 conv(1 1 1, 1 1 1) 执行后即得到y(n)为 ans =1 2 3 2 1 注意：使用conv 函数时，h(n) 和x(n)都必须是有限长的, ，否则不能使用conv 函数。 例5-1 时域离散序列的卷积计算 与图形显示 例5-1 (教材p63) ：已知离散信号x(n)和h(n)， 求y(n)=x(n)*h(n)，并用图形表示。例5-1的MATLAB程序Nh=20;Nx=10;m=5;%设定Nx，Nh和位移值m n=0:Nh-1;h1=(0.9).n;%产生h1(n) h2=h1; nx=0:Nx-1;x1=ones(1,Nx);%产生x1(n) x2=zeros(1,Nx+m); for k=m+1:m+Nx%产生x2(n)=x1(n-m) x2(k)=x1(k-m); end %产生x2(n) y1=conv(x1,h1);%计算y1(n)=x1(n)*h1(n) y2=conv(x2,h2);%计算y2(n)=x2(n)*h2(n) subplot(3,2,1) stem(nx,x1,'.') axis(0 30 0 1.2),title(x1(n) %绘图 （以下省略）5.2.2 滤波 数字滤波器的系统函数H(z)用如下式表示: 在MATLAB中，用向量b，a来表示滤波器的系 数b(i)和 a(i)。 滤波器分类 当n = 0，m0时，称为AR滤波器，即自回归（ Auto Recurrence）滤波器，具无限冲激响应（ IIR），也即其单位采样响应h(n)具无限长度； 若m = 0，a(1) 0，称为MA滤波器，即滑动平 均（Moving Average）滤波器，其单位采样响 应h(n)是有限长度，故称有限冲激响应（FIR） 滤波器； 如果n、m都大于零，称为ARMA滤波器，而其 冲激响应也为IIR。filter函数 MATLAB提供了 filter函数来对离散信号进行 滤波，表达信号通过系统后的结果。 与conv不同的是，filter函数可适用于无限冲 激响应系统的情况，但信号仍须是有限长的 。 例如，一个单极点的低通滤波器系数如下： >>b = 1; % 分子系数向量b(i)>>a = 1 -0.9; % 分母系数向量a(i)如果用filter函数实现对信号x滤波，只要调用 ： >>y = filter(b,a,x);就可给出输入x经过滤波以后的输出y。5.2.3 单位冲激响应 数字滤波器的单位冲激响应定义为输入 为单位样本序列时数字滤波器的响应 ， 即： h(n) = T (n) 其中：单位冲激响应的MATLAB实现 MATLAB近似实现单位采样信号的方法为： imp = 1; zeros(p,1); % zeros(p,1)产生 p个零元素组成的列向量，p是 正整数。使用imp后，滤波器的冲激响应可近 似得到为： h = filter(b,a,imp); impz 函数可以直接求出数字滤波器的单位冲激 响应，即：impz(b,a) 该命令将同时绘出滤波器的单位冲激响应 ，教材p66图5-2。53 频域和Z域分析 频率响应 ，零极点分析 5.3.1 频率响应 MATLAB数字信号处理工具箱有很多函数提供 对模拟和数字滤波器的频率响应分析。其中， freqz 函数和freqs 函数分别返回数字和模拟滤 波器的频率响应。 工具箱中通常使用的单位频率是Nyquist频率， 即采样频率的1/2。 注意：就数字滤波器函数来说，其频域指标中 的所有频率都以Nyquist频率进行归一化。 因此Nyquist频率也称归一化频率。 关于Nyquist频率的说明 例如：系统采样频率为1000Hz，则若数字滤波 器的截止频率等于300Hz，经Nyquist频率归一 化后，其归一化频率就是300/500=0.6。若将归 一化频率转换成数字信号处理教科书中所使用 的数字频率(rad)，需乘以；反之，若乘以采 样频率fs的一半，则将归一化频率转换回了模 拟域频率(Hz)。 归一化频率f 应满足0<f<1。 对于频率响应而言，归一化频率和模拟频率都 可使用。freqz命令 freqz 函数提供了基于FFT算法所得到的数字滤波器 H(z)的频率响应 最常用的形式： h,w = freqz(b,a,p) 其意义是，对于数字滤波器：freqz 函数接受滤波器系数向量b和a，以及一个用于指 定所需计算的频率响应样点数的整数p，返回一个与由 p个频率样本点构成的实频率向量相对应的复频率响应 向量h。 freqz的命令形式 h = freqz(b,a,w) 采用上面的形式时，需先对频率样本点向量w作出 定义。通常的做法是使用linspace函数。 h,w = freqz(b,a) w和p没有定义 ，默认w由（0）上均分的512点构 成，频率单位为rad/sample。 h,w = freqz(b,a,p，whole) 使用带参数whole选项的命令形式 h,w = freqz(b,a,p,fs) 用参数选项指定采样频率fs linspace的命令形式 linspace函数的命令形式有两种： linspace（x1,x2） %产生x1,x2范围内经线性分割得到的100点的向 量 linspace（x1,x2，N） %产生x1,x2范围内经线性分割得到的N点的向量 例如： w = linspace（0，pi） % w由（0，pi）上100个等分点组成w = linspace（0，1000） % w由（0，1000）上100个等分点组成freqz函数无返回输出参数格式的 调用 freqz函数还能以无返回输出参数的格式调用： freqz（b，a，256）或 freqz（b，a，256，2000） 无返回输出参数调用freqz函数的好处是，MATLAB能 够自动绘出频率在（0）范围内的幅频特性和相频特 性图。 注意，无返回输出参数调用freqz函数绘出的相频特性 不能正确给定在 处的值，这是因为在MATLAB中， 属于下半个单位圆。 频率响应的实例 例：先构成一个截止频率为400Hz的9阶巴特沃思（ Butterworth）低通数字滤波器，求出其系数b，a，再求 出其256点频率响应。指定的采样频率fs =2000Hz。 实现1：先调用butter函数，再调用freqz函数； 实现2：无返回输出参数，调用freqz函数； 实现3：为得到完整的（0）范围内的幅相特性图， 可用带返回输出参数的格式调用freqz，在得到各频率 点上的频率响应的数值后，使用MATLAB所提供的abs 函数和angle函数分别提取幅频特性与相频特性。 注意：为了得到连续的相频曲线，需要使用unwrap函数，该 函数通过在发生跳变后的各处自动加上，从而除去了相位的跳变，这称为相位的“解卷绕”。 实现1b,a = butter(9,400/1000); % 括号中第一个输入是阶数，第二个是归一化 截止频率 h,f = freqz(b,a,256,2000);或 h = freqz（b，a，256，2000）；实现2 实现1中第二条命令的形式 可改写为freqz（b，a，256） 或freqz（b，a，256，2000 ） 可自动绘出频率在（0 ）范围内的幅频特性和相 频特性图。 注意无返回输出参数调用 freqz函数绘出的相频特性 不能正确给定在=处的 值，因为= 属于下半个 单位圆。 实现3b,a = butter(9,400/1000); w = linspace(0,1000); h = freqz(b,a,w,2000); mag = abs(h); pha = angle(h); % 提取滤波器频率响应的幅度mag和相位pha semilogy(w,mag); title('Magnitude'); figure; plot(w,pha); title('Phase'); 出现了幅度 为2的跳变 卷绕解卷绕： unwrap函数 为了得到连续的相频曲线，需要使用unwrap函数，该 函数通过在发生2跳变后的各处自动加上2 ，从而 除去了相位的跳变，这称为相位的“解卷绕”。解卷绕freqs函数 freqs函数用于求取由b和a系数向量定义的模拟滤 波器 H（s）的频率响应。 使用方法类似freqz函数。 与第二章（p32例2-1）采用数组相除方法求取频 率响应相比，使用freqs 函数要方便很多。 5.3.2 零极点分析 zplane 函数用于画出线性系统在Z平面上的零 极点。有两种使用方法： 1、在已知零极点时，例如某滤波器的零点为- 1/2，一对共轭极点为 和 时， 只要输入命令 zer = -0.5; pol = 0.9*exp(j*2*pi*-0.3 0.3'); zplane(zer,pol) 即可画出零极点。 （见p70图5-6）用zplane函数绘制零极点图 另一种情况：已知系统的系统函数系数向量b 和 a ，则可通过调用zplane(b,a)绘出零极点。 这种情形下，zplane 函数先求得系统函数的零 点和极点，然后绘出零极点图。5.4 滤波器设计滤波器技术指标， IIR滤波器设计，FIR滤波器设计5.4.1 滤波器的技术指标 设计滤波器之前，要定义一组滤波器的技术指标。 滤波器的技术指标通常采用幅度指标给出，以保证物 理上的可实现性。 幅度指标以容差方式给出 ，如p71图5-7所示。 具体指标：通带截止频率Wp ；阻带截止频率Ws ；通 带波纹p ；阻带衰减s。其中通带波纹和阻带衰减 也可以用分贝数给出 滤波器以容差方式给出的幅度指 标 |H( jw)|10通带过渡带阻带pwswsdpd-1w5.4.2 IIR 滤波器设计 MATLAB工具箱提供了几种模拟滤波器原型的 产生函数,如巴特沃思(Butterworth),切比雪夫 (Chebyshev)滤波器等。 还提供了从模拟低通滤波器原型转换为高通、 带通和带阻的转换函数，模拟滤波器转换为数 字滤波器的双线性变换法和冲激响应不变法， 模拟IIR滤波