数字信号处理实验考试MATLAB笔记

数字信号处理实验考试 MATLAB 笔记 1 计算卷积 y conv x h 2 filter 函数计算差分方程的零状态相应 离散 LTI 都可用如下的线性常系数的差分方程描述 N n M n nn nkxbnkya 00 filter 函数 可方便地计算出上述差分方程的零状态响应 filter 函数调用形式为 y filter b a x a 对应 y b 对应 x x 表示输入的序列 其中 a a0 a1 aN b bo b1 bM h n impz b a L 是求单位脉冲响应的函数 L 是点数 3 freqz 函数计算有理分式的频率相应的抽样值 若 X 可表示为 j e j eX j j eA eB Njj Mjj aNeeaa bMeebbo 10 1 freqz 的调用形式为 X freqz b a w 其中 b b0 b1 bM a a0 a1 aN W 为抽样的频率点 abs angle real imag 分别计算 X 的幅度 相角 实部 和虚部 4 X z 的部分分式展开的函数 它的调用形式如下 r p k residuez b a 分解为并联分解为并联 1 型 型 residue 分解为并联分解为并联 2 型 型 要弄清楚它们的用法要弄清楚它们的用法 b a 分别表示用 z表示 X z 的分子和分母多项式 注意 b b0 b1 b1 bn 1 1 其中 bn 是 z n 的系数 如果 X z 的部分分式展开为 X z 1 21111 21 31 4 31 3 21 2 11 1 zkk zp r zp r zp r zp r 则 residuez 的返回参数 r p k 分别为 r r r r r 1234 p p1 p2 p3 p3 k k1 k2 同一极点 p3 在向量 p 中出现了两次 表示 p3 是一个二阶的重极点 Residuez 也用于由 r p k 计算 z表示的 X z 的分子和分母多项式 其调用形式为 1 b a residuez r p k 5 tf2zp 和和 zp2tf 对系统函数 H z 的形式转换 z 正幂有理多项式表示的系统函数为 2 1 2 1 1 1 Nazaza Mbzbzb zH NN MM 用零点 极点和常数表示的一阶因子形式的系统函数为 2 1 2 1 Npzpzpz Mzzzzzz kzH 其中 z p k tf2zp b a 和 b a zp2tf z p k 其中 a a 1 a N b b 1 b M 注意缺失项要补零 zplane b a 用于给出参数直接画图 画出 z 平面的零极点分布来判断系统 的稳定性 6 利用 利用 MATLAB 计算信号计算信号 DFT 在 MATLAB 信号处理工具箱中 函数 dftmtx N 可用来产生 N N 的 DFT 矩正 D N N 的 IDFT 矩正 D可用函数 conj dfmtx N N 来确定 此外 N 1 N MATLAB 提供了 4 个内部函数用于计算 DFT 和 IDFT 它们分别是 fft x fft x N ifft X ifft X N fft x 计算 M 点的 DFT M 是序列 x 的长度 即 M length x 先对 x 离散 才能 DFT 变换 fft x N 计算 N 点的 DFT 若 M N 则将原序列截短为 N 点序列 再计 算其 N 点 DFT 若 MN 则将原序列截短为 N 点序列 再计 算其 N 点 IDFT 若 M N 则将原序列补零至 N 点 然后计算其 N 点 IDFT 7 用 用 MATLAB 实现模拟低通滤波器的设计实现模拟低通滤波器的设计 MATLAB 信号处理工具箱提供了常用的设计模拟低通滤波器的 MATLAB 函数 在实现应用中 可方便地调用这些函数完成模拟滤波器的设计 关于这 些函数现分别介绍如下 Butterworth 滤波器滤波器 n wc buttord wp ws rp rs s num den butter n wc s 下面这个形式是设计下面这个形式是设计IIR数字滤波器 其中数字滤波器 其中wp 2 fp ft ws 2 fs ft 与模拟滤波器与模拟滤波器 设计形式不同 设计形式不同 wn wp ws都是归一化的数值都是归一化的数值 N Wn buttord Wp Ws Rp Rs b a butter N Wn 根据 BW 型滤波器的设计指标 利用 MATLAB 函数 buttord 获得 BW 型 滤波器参数 N 和 wc 函数 buttord 的输入参数 wp 和 ws rad s 分别表示滤波器 的通带和阻带截频 Ap 和 As dB 表示滤波器的通带和阻带衰减 s 表示所设 计的是模拟滤波器 函数 buttord 返回参数 N 为 BW 滤波器的阶数 wc rad s 等 于 BW 滤波器 3 dB 截频 由于 wc 由阻带方程确定 故由参数 N wc 得出 c 的滤波器在阻带刚好满足设计指标 在通带将超过设计指标 当 BW 型滤波器的参数 N wc 确定后 可用 MATLAB 函数 butter 获得 BW 型滤波器系统函数的分子多项式 num 和分母多项式 den sHLP Chebyshev I 型滤波器型滤波器 n wc cheb1ord wp ws rp rs s num den cheby1 n rp wc s MATLAB 函数 cheblord 返回参数 N 表示 CB I 型滤波器的阶数 wc rad s 等于 wp 参数 N wc 的取值可使 cheby1 设计出的 CB I 型滤波器在通带刚好 满足设计指标 cheby1 函数利用参数 N wc 和 Ap 确定 CB I 型滤波器系统函数的 sHLP 分子多项式 num 和分母多项式 den Chebyshev II 型滤波器型滤波器 N wc cheb2ord wp ws Ap As s num den cheby2 N As wc s MATLAB 函数 cheb2ord 返回参数 N 表示 CBII 型滤波器的阶数 wc rad s 的取值可使 cheby2 设计出的滤波器在通带刚好满足设计指标 cheby2 函数利用参数 N wc 和 As 确定 CBII 滤波器系统函数的分 sHLP 子多项式 num 和分母多项式 den 椭圆滤波器椭圆滤波器 N wc ellipord wp ws Ap As s num den ellip N Ap As wc s 数字滤波器设计的情况 n wc ellipord wp ws ap as num den ellip n ap as wc MATLAB 函数 ellipord 返回参数 N 表示椭圆滤波器的阶数 wc wp MATLAB 函数 ellip 确定阶数为 N 通带衰减为 Ap dB 阻带衰减为 As dB 通带截频为 wc 的椭圆滤波器系统函数的分子多项式和分母多项式 在 sHLP 滤波器的阶数确定后 采用了 4 1 3 节中方案 A 重新计算滤波器的参数 k1 故 设计出的椭圆滤波器阻带衰减超过设计指标 而其他指标刚好满足设计要求 8 用用 MATLAB 实现模拟域频率变换实现模拟域频率变换 MATLAB 信号处理工具箱提供了实现四种模拟域频率变换的函数 它们 分为 低通到低通的变换 numt dent 1p21p num den w0 低通到高通的变换 numt dent 1p2hp num den w0 低通到带通的变换 numt dent 1p2bp num den w0 B 低通到带阻的变换 numt dent 1p2bs num den w0 B 其中 num den 分别表示变换前模拟滤波器系统函数的分子多项式和分母多项 式 numt dent 分别表示变换后模拟滤波器系统函数的分子多项式和分母多项 式 w0 和 B 为变换中的参数 9 脉冲响应不变的脉冲响应不变的 MATLAB 实现实现 MATLAB 信号处理提供的 impinvar num den Fs 函数 可实现用脉冲响应 不变法将模拟滤波器转化为数字滤波器 起调用形式为 numd dend impinvar num den Fs 式中 num 和 den 分别表示模拟滤波器系数函数 H s 的分子多项式和分母多项 式 Fs 是脉冲响应不变法中的抽样频率 单位是 Hz 输出变量 numd 和 dend 分别表示数字滤波器的系统函数 H z 分子多项式和分母多项式 若某个因果模拟滤波器的系统函数为 52 142 2 ss s sH 取 T 1 则由 numd dend impinvar 2 14 1 2 5 1 得 numd 2 0000 2 3133 dend 1 0000 0 3062 0 1353 所以由脉冲响应不变法获得的数字滤波器为 1353 0 3062 0 1 3133 2 2 1 1 z z zH 10 双线性变换法的双线性变换法的 MATLAB 实现实现 MATLAB 信号处理工具箱提供的 num den Fs 函数可以用来实现双线性 变换 其调用形式为 numd dend bibinear num den Fs 式中 num 和 den 分别表示模拟滤波器系统函数 H s 的分子多项式和分母多项 式 Fs 1 T 输出变量 numd 和 dend 分别数字滤波器系统函数 H z 的分子多项式和 分母多项式 例例 4 19 三阶归一化 Butterworth 滤波器的系统函数为 1 1 1 2 sss sH 解 取 T 2 则由 den conv 1 1 1 1 1 numd dend bilinear 1 den 0 5 得经双线性变换后的数字滤波器的分子多项式和分母多项式为 numd 0 1667 0 5000 0 5000 0 1667 dend 1 0000 0 0000 0 3333 0 0000 注意 0 1667 是 1 6 的近似值 0 3333 是 1 3 的近似值 所以双线性变换后 的数字滤波器的系统函数 H z 为 2 31 3 1 1 1 6 1 z z zH 11 用 用 MATLAB 实现数字滤波器设计实现数字滤波器设计 MATLAB 也提供了直接设计 IIR 数字滤波器的函数 设计的基本思想是 用式 4 103 将数字滤波器的频率指标转化为模拟滤波器的指标 然后设计模 拟滤波器 最后用双线性变换把模拟滤波器转换为数字滤波器 由于在设计中 用的模拟滤波器的类型不同 所以给出了以下四组不同的函数 BW 型数字低通滤波器型数字低通滤波器 N Wc buttord Wp Ws Ap As num den butter N Wc 在函数 buttord 函数中 调用参数 Wp Ws 是数字低通滤波器的归一化的 通带和阻带截频 例如要求数字滤波器的 则 4 0 1 0 sp Wp 0 1 Ws 0 4 Ap 和 As 为滤波器的通带和阻带的衰减 dB 返回的参数为 和 Wc 其中 N 为滤波器阶数 函数 butter 中 调用参数 N 和 Wc 由函数 buttord 确定 返回参数 num 和 den 分别是数字滤波器的分子多项式和分母多项式的系数 CB I 型数字低通滤波器型数字低通滤波器 N Wc cheb1ord Wp Ws Ap As num den cheby1 N Ap Wc CB II 型数字低通滤波器型数字低通滤波器 N Wc cheb1ord Wp Ws Ap As num den cheby1 N Ap Wc 椭圆型数字低通滤波器椭圆型数字低通滤波器 N Wc ellipord Wp Ws Ap As num den ellip N Ap As Wc 利用上述四组 MATLAB 函数 也可设计数字高通 带通和带阻滤波器 调用方法可参见 MATLAB 手册 12 利用 利用 MATLAB 实现实现 FIR 滤波器设计滤波器设计 5 5 1 窗函数法的窗函数法的 MATLAB 实现实现 窗函数的计算窗函数的计算 MATLAB 提供了许多常用的窗函数 其中部分窗函数的调用形式为 w hanning N w hamming N w blackman N w Kaiser N beta 其中 N 是窗函数的长度 beta 是控制 kaiser 窗形状的参数 返回的变量 w 是一个 长度为 N 的列向量 给出窗函数在 N 点的取值 窗函数法设计窗函数法设计 FIR 滤波器滤波器 窗函数法设计 FIR 滤波器一般分为 3 个步骤 第 1 步估计 FIR 滤波器阶数 M 或长度 N 如果用 Kaiser 窗时可用式 5 41 估计 FIR 滤波器阶数 第 2 步确定所用的窗函数并计算出窗函数的值 第 3 步计算理想滤波器的单位脉冲 响应并用窗函数将其截断即得所设计的 FIR 滤波器的 h k 例例 5 5 用 I 型线性相位滤波器设计一满足下列指标的 FIR 高通滤波器 dBAsdBApradrad sp 40 3 0 7 0 8 0 解 表 5 2 可知用 Hann 窗可使滤波器满足指标 由过渡带宽宽度得滤波 器长度 N 需满足 62 2 6 sp N 由于要求用 I 型滤波器 取 N 63 取理想高通滤波器的截频为 radspc 75 0 2 由式 5 28 得 31 75 0 75 0 31 kSkkh ad 用长度 N 63 的 Hamming 窗截断即得所要求的 FIR 高通滤波器的 h k 为 khd kwkhkh d 用 Kaiser 窗得出的滤波器阶数较低 用 Hann 设计出的滤波器阻带的波纹 衰减较快 用 Kaiser 窗设计 FIR 滤波器阶数的 MATLAB 函数基本调用形式为 M Wc beta ftype Kaiserord f a dev h fir1 M Wc ftype window 函数 Kaiserord 完成估计滤波器阶数 M 及 Kaiser 窗的参数 beta 函数 fir1 根据 所用窗函数及理想滤波器的截频计算出 FIR 滤波器的 h 函数 Kaiserord 调用参数 f 表示需设计的 FIR 滤波器的频带 如 FIR 滤波 器的 B 个频带分别为 221 1 22 1 254 321 0 BBB fffff fff 则 f 是一个有 2B 2 个元素的向量 B2 其值为 2254321 B fffffff 函数 kaiserord 调用参数 a 是一个 B 个元素的向量 分别表示 FIR 滤波器 在B 个频带中的幅度值 一般对通带取值为 1 阻带取值为 0 函数 kaiserord 调用参数 dev 是一个 B 个元素的向量 分别表示 FIR 滤波 器在 B 个频带中的波动值 由前面的分析可知 Kaiser 窗设计的 FIR 滤波器在 各带的波动是相等的 所以 kaiserord 是利用最小的波动值来估计滤波器的阶数 函数 kaiserord 返回参数 M 及 beta 分别表示 FIR 滤波器阶数 M 及 kaiser 窗 的参数 beta 返回参数 Wc 和 ftype 是函数 firl 的调用参数 Wc 是 B 1 个元素 的向量 如 Wc W1 W2 W 1 B 则表示理想 FIR 滤波器的 B 个频带分别为 132211 0 B wwwwww 返回参数 ftype 是一个字符串 ftype 为空时 表示滤波器为低通 ftype high 滤波器为高通 ftype stop 滤波器为带阻 ftype DC 0 多带滤波器第一个 频带为阻带 ftype DC 1 多带滤波器第一个频带为通带 例如 f 0 2 0 6 a 1 0 时 函数 kaiserord 返回值 Wc 和 ftype 分别为 Wc 0 4 ftype 当 f 0 2 0 6 a 0 1 时 函数 kaiserord 返回值 Wc 和 ftype high 当 f 0 2 0 4 0 6 0 8 a 1 0 1 时 函数 kaiserord 返回值 Wc 和 ftype 分别为 Wc 0 4 ftype stop 当 f 0 2 0 4 0 6 0 8 a 0 1 0 时 函数 kaiserord 返回值 Wc 和 ftype 分别为 Wc 0 3 0 7 ftype DC 0 函数 firl 的调用参数 M 表示滤波器的阶数 Wc 表示理想 FIR 滤波器的 B 个频带 ftype 表示滤波器的类型 缺省值为空白 window 是一个长度为 M 1 的向量 如果调用是没给窗函数 程序 firl 自动使用 Hamming 窗 MATLAB 信号处理工具箱提供的另一个用窗函数法设计 FIR 滤波器的函数 为 h fir2 M f a npt window fir2 实现的基本思想是利用线性内插得出理想滤波器的频率响应 然后 j eH 对均匀抽样得出在 npt 个点 缺点值为 512 上的抽样值 对 j eH j eH 的 npt 点的 IDFT 得 h k 最后用窗函数将 h k 截断为长度为 M 1 的 j eHdd 序列 函数 fir2 的调用参数 N 表示的阶数 npt 表示 IDFT 的点数 缺省值为 512 window 为所用窗函数序列 缺省时使用 Hamming 窗 如果所设计的滤波 器的 B 频率点上的幅度值为 1 0 21 2 Bm fmj fffaeH 则 f 和 a 的值分别为 f f1 f2 f a a1 a2 a BB 函数 fir2 的返回参数 h 为长度为 M 1 的 FIR 滤波器系数 13 数字滤波器结构的数字滤波器结构的 MATLAB 实现实现 7 4 1 直接型直接型 IIR 直接型结构由两个矢量 b 和 a 描述 b 包含 bn 系数 a 包含 an 系数 它由 filter 函数实现 调用格式为 filter b a x FIR 直接型结构由包含 bn 系数的矢量 b 描述 它也由 filter 函数实现 只 是把矢量 a 置为 1 既调用格式为 filter b 1 x 线性相位 FIR 数字滤波器的数字滤波器的系统函数形式相同 只是系数具 有对称关系 因此 在 MATLAB 实现上 线性相位结构可以由直接型结构实 现 7 4 2 级联型级联型 对于 IIR 级联型实现 系统函数 H z 通常可表示为有限个实系数二阶有理 分式之积 即 7 41 2 2 1 10 2 2 1 10 1 zazaa zbzbb zH iii iii L i MATLAB 用一个 L 6 矩阵 sos 来表示这个二阶分式 sos 7 42 LLLLLL aaabbb aaabbb aaabbb 221210 221202221202 121101121101 若系统函数 H z 是以零极点增益表示的 即 7 43 21 21 n n pzpzpz zzzzzz kzH 则可用 zp2cos 函数将其转化式 7 41 表示的二阶分式之积 调用格式为 sos zp2sos z p k 将式 7 41 中的分母设为 1 即得 FIR 系统函数的有限个二阶因式乘积表示 因此 FIR 系统函数的每一个二阶因式也可用 sos 来表示 只需设式 7 42 中 的 a 1 其余 a 系数为零即可 i0 例例 7 5 实现下面系统的级联型结构 321 21 6 1 3 1 6 1 1 3 2 3 5 3 zzz zz zH 解 program 7 1 Cascade Realization of an IIR format rat num 3 5 3 2 3 den 1 1 6 1 3 1 6 z p k tf2zp num den sos zp2sos z p k 程序运行结果为 sos 2 3 0 0 1 3 0 9 2 5 2 1 1 1 2 1 2 由此可得实系数有理分式之积的形式 21 21 1 2 1 2 1 1 2 5 2 9 3 1 1 3 2 zzz zz z zH 一般 系统函数的实数极点对应系数一阶有理分式 共轭复数极点对应系数二 阶有理分式 7 4 3 并联型并联型 对于 IIR 并联型实现 系统函数 H z 通常可以表示为有限个实系数一阶或 二阶有理分式之和的形式 这种形式在 MATLAB 中可以通过 residuez 函数实现 如例 7 5 所述系统 可以下面程序将 H z 展开成部分分式之和 program 7 2 parallel Realization of an IIR format rat num 3 5 3 2 3 den conv 1 1 3 1 1 2 1 2 r p k residuez num den 程序运行结果为 r 1 2 1012 1785i 1 2 1012 1785i 2 p 1 4 506 765i 1 4 506 765i 1 3 k 该系统有一对共轭复数极点 利用 b1 a1 residuez R1 P1 0 语句可获得其所对 应的实数二阶分式的分子 分母多项式系数 其中 R1 为共轭复数留数所构成 的向量 P1 为共轭复数极点所构成的向量 b1 a1 为有理分式分子和分母多项 式的系数向量 运行下面语句 R1 r 1 r 2 P1 p 1 p 2 b1 a1 residuez R1 P1 0 可得实系数二阶因子式的分子 分母多项式系数分别为 b1 1 1 0 a1 1 1 2