1、第7章 多采样率数字信号处理,7.1 引言 7.2 学习要点及重要公式 7.3 采样率转换系统的高效实现 7.4 教材第8章习题与上机题解答,7.1 引 言本章内容与教材第8章内容相对应。在实际系统中， 经常会遇到采样率的转换问题， 要求一个数字系统能工作在“多采样率”状态。 这样的系统中， 不同处理阶段或不同单元的采样频率可能不同。 例如， 在DSP开发仿真实验系统中， 为了抗混叠滤波器设计实现简单， 降低系统复杂度， 应先统一对模拟信号以系统最高采样频率采样， 然后， 根据实验者选择的各种采样频率， 在数字域改变采样频率。,由于在数字域改变采样频率完全用软件实现， 所以， 采样频率可以任意
2、选择， 又避免了对各种不同的采样频率设计各种不同的抗混叠模拟滤波器， 从而使系统复杂度大大降低。 采样率转换在现代通信、 信号处理和图像处理等领域都有广泛的应用， 应用实例举不胜举， 本科大学生能理解的几种实例在教材8.1节有介绍， 这里不再重复。,7.2 学习要点及重要公式本章要求学生熟悉采样率转换的基本概念和种类， 了解采样率转换的应用价值和适用场合， 掌握三种常用的采样率转换基本系统（整数因子D抽取、 整数因子I插值和有理数因子I/D采样率转换）的基本原理、 构成原理方框图及其各种高效实现方法（FIR直接实现、 多相滤波器实现)， 每种实现方法的特点。 这些专业基础知识对进一步学习、 设
3、计、 开发工作在多采样率状态的各种复杂系统是极其重要的。 采样率转换的基本概念和种类， 以及应用价值和适用场合在教材中已有较详细的介绍， 这里不再重复。,下面对三种常用的采样率转换基本系统的重要知识点及相关公式进行归纳总结， 以便读者复习巩固。 值得注意， 要理解采样率转换原理， 必须熟悉时域采样概念、 时域采样信号的频谱结构、 时域采样定理。 此外， 时域离散线性时不变系统的时域分析和变换（Z变换、 傅里叶变换）域分析理论是本章的分析工具。 只有熟练掌握上述基础知识， 才能掌握本章的知识， 否则， 无法理解本章内容。,7.2.1 整数因子D抽取1. 整数因子D抽取器原理框图按整数因子D对x(
4、n)抽取的原理方框图如图7.2.1所示。,图7.2.1,2. 整数因子抽取器的功能整数因子抽取器的功能表现为： 输出端信号采样频率Fy降为输入端信号采样频率Fx的1/D， 即Fy=Fx/D。 3. 知识要点及重要公式经过抽取使采样率降低， 会引起新的频谱混叠失真， 所以， 必须在抽取前进行抗混叠滤波。 抗混叠滤波器hD(n)有两项指标： 根据采样定理， 抗混叠滤波器的阻带截止频率为,相应的数字阻带截止频率为,(7.2.1),(7.2.2),应当注意， 由于抗混叠滤波器工作于输入信号采样频率Fx， 所以， 式（7.2.2）中用Fx换算得到相应的数字截止频率， 绝对不能用Fy换算。,抗混叠滤波器的
5、通带截止频率（或过渡带宽度）取决于抽取系统对信号频谱的失真度要求。 设计时根据具体要求确定抗混叠滤波器的其他三个指标参数（通带截止频率p、 通带最大衰减p、 阻带最小衰减s）。 例如， 要求抽取过程中频带0, fp上幅频失真小于1%（显然fpfs）， 由于抽取引起的频谱混叠失真不超过0.1%。 这时， 抗混叠滤波器的通带截止频率为p=2fp/Fx， 通带最大衰减为p=20 lg(11%)=0.0873 dB， 阻带最小衰减为s=20 lg(0.1%)=60 dB。,整数因子抽取器的输入x(n)和输出y(m)的关系式如下：,(7.2.3),(7.2.4),(7.2.5),在主值区, 上Y(ejy
6、)为,|y|,7.2.2 整数因子I内插器1. 整数因子I内插器原理框图按整数因子I对x(n)内插的原理方框图如图7.2.2所示。,图7.2.2,2 整数因子内插器的功能整数因子内插器的功能表现为： 输出端信号采样频率Fy提高为输入端信号采样频率Fx的I倍， 即Fy=IFx。 3. 知识要点及重要公式经过整数倍零值内插使采样率升高至I倍， 不会引起新的频谱混叠失真， 但是会产生I-1个镜像频谱（见图7.2.3(b)）。 从时域考虑， 零值内插器输出v(m)的两个非零值之间有I-1个零样值， 不是我们所期望的提高采样率后的采样序列y(m)。 所以， 必须进行低通滤波， 滤除镜像频谱， 得到y(m
7、)。,镜像频谱滤波器hI(n)的指标也有两项： 阻带截止频率和数字域阻带截止频率。 根据采样理论， 采样频率提高到Fy时， 采样信号序列y(m)的频谱以Fy为周期。 由于输入端信号x(n)的采样频率为Fx， 所以x(n)的频带宽度不会超过Fx/2（对应的数字频率为）。 因此， 整数因子I内插器输出信号y(n)的频谱一定是带宽为Fx/2、 重复周期为Fy的周期谱。 我们将零值内插器输出的信号v(m)的频谱中除了y(m)频谱以外的其他频谱称为镜像频谱。 镜像频谱滤波器的作用就是让y(m)频谱尽量无失真通过， 滤除v(m)中的镜像频谱。 当I=3时， 图7.2.2中各点信号的频谱示意图如图7.2.3
8、所示。,图7.2.3,综上所述可知， 镜像频谱滤波器的阻带截止频率为,相应的数字域阻带截止频率为,(7.2.6),(7.2.7),应当注意， 由于镜像频谱滤波器工作于输出信号采样频率Fy， 所以， 式（7.2.7）中用Fy换算得到相应的数字截止频率， 绝对不能用Fx换算。,与抗混叠滤波器情况类似， 根据具体要求确定镜像 频谱滤波器其他三个指标参数（通带截止频率p、 通带最大衰减p、 阻带最小衰减s）。 镜像频谱滤波器的p（或过渡带宽度）和p取决于内插系统对有用信号频谱的保真度要求， 而s取决于内插系统对镜像高频干扰的限制指标。理想情况下， 镜像频谱滤波器的频率响应函数为,(7.2.8),通带内
9、幅度常数取I是为了确保在m=0, I, 2I, 3I, 时， 输出序列y(m)=x(m/I)。 整数因子I内插系统的时域输入输出关系式如下：,(7.2.9),因为除了在I的整数倍点v(kI)=x(k)以外， v(k)=0， 所以,(7.2.10),变换域输入输出关系式：,(7.2.11),计算单位圆上的V(z)得到v(m)的频谱为,(7.2.12),(7.2.13),(7.2.14),理想情况下， HI(ej)由式（7.2.8）确定， 所以，,(7.2.15),图7.2.4,7.2.3 有理数因子I/D采样率转换系统1 有理数因子I/D采样率转换系统原理框图有理数因子I/D采样率转换的原理框图
10、如图7.2.4所示。,2 有理数因子I/D采样率转换系统的功能有理数因子I/D采样率转换系统首先对输入序列x(n)按整数因子I内插， 然后再对内插器的输出序列按整数因子D抽取， 达到按有理数因子I/D的采样率转换。 如果仍用Fx=1/Tx和Fy=1/Ty分别表示输入序列x(n)和输出序列y(m)的采样频率， 则Fy=(I/D)Fx。应当注意， 先内插后抽取才能最大限度地保留输入序列的频谱成分(请读者解释为什么)。,3 知识要点及重要公式图7.2.4中滤波器h(l)同时完成镜像滤波和抗混叠滤波功能。 所以， 理想情况下， 滤波器h(l)是理想低通滤波器， 其频率响应为,图7.2.4中各点信号的时
11、域表示式归纳如下：,(7.2.16),线性滤波器输出序列为,(7.2.17),线性滤波器输出序列为,(7.2.18),整数因子D抽取器输出序列y(m)为,(7.2.19),式（7.2.19）就是有理数因子I/D采样率转换系统的输入输出时域关系。 如果线性滤波器用FIR滤波器实现， 则式（7.2.19）为有限项之和， 所以可以直接按式（7.2.19）编程序计算输出序列y(m)。 当然， 也可以采用教材上介绍的各种高效实现结构以硬件或软硬结合来实现。,7.3 采样率转换系统的高效实现 实际上， 采样率转换系统的高效实现就是指其中的FIR数字滤波器的高效实现。 这里高效的含义有三个方面： 在满足滤波
12、指标要求的同时， 滤波器的总长度最小； 使滤波处理计算复杂度最低； 对滤波器的处理速度要 求最低。 教材中介绍了采样率转换系统的两种实现方法： 直接型FIR滤波器结构、 多相滤波器实现。 各种实现方法的原理、 结构及其特点在教材中都有较详细的叙述， 本书不再重复。,7.4 教材第8章习题与上机题解答1. 已知信号x(n)=anu(n), |a|1。 （1） 求信号x(n)的频谱函数X(ej)=FTx(n)； （2） 按因子D=2对x(n)抽取得到y(m)， 试求y(m)的频谱函数。 （3） 证明： y(m)的频谱函数就是x(2n) 的频谱函数。 解：,（2） 根据式（7.2.5a）可知y(m)
13、的频谱函数为,（3）,2. 假设信号x(n)及其频谱X(ej)如题2图所示。 按因子D=2直接对x(n)抽取， 得到信号y(m)=x(2m)。 画出y(m)的频谱函数曲线， 说明抽取过程中是否丢失了信息。,题2图,解： 由抽取原理知道， 按因子D=2对x(n)抽取时， 抗混叠滤波器的截止频率为/D rad， 所以得到信号y(m)=x(2m)丢失原信号x(n)中|/D的频率成分信息。 本题中，D=2， 而x(n)的带宽为/3， 所以， 抽取过程中不会丢失信息。 由时域采样理论， 画出y(m)的频谱如题2解图所示。,题2解图,3. 按整数因子D=4抽取器原理方框图如题3图（a）所示。 其中， Fx
14、=1 kHz, Fy=250 Hz, 输入序列x(n)的频谱如题3图（b）所示。 请画出题3图（a）中理想低通滤波器 hD(n)的频率响应特性曲线和序列v(n)、 y(m)的频谱特性曲线。,题3图,解： 抽取因子D=Fx/Fy=1000/250=4。 题3图（a）中理想低通滤波器hD(n)的频率响应特性曲线以及序列v(n)和y(m)的频谱特性曲线分别如题3解图（a）以及（b）和（c）所示。,题3解图,4. 按整数因子I内插器原理方框图如题4图所示。 图中， Fx=200 Hz, Fy=1 kHz， 输入序列x(n)的频谱如题3图（b）所示。 确定内插因子I， 并画出题4图中理想低通滤波器hI(
15、n)的幅频响应特性曲线和序列v(m)、 y(m)的频谱特性曲线。,题4图,解： 内插因子I=Fy/Fx=1000/200=5。 题4图中序列v(m)、 理想低通滤波器hI(n)的幅频响应特性曲线和y(m)的频谱特性曲线分别如题4解图（a）、 （b）和（c）所示。 5*. 设计一个抽取器， 要求抽取因子D=5。 用remez函数设计抗混叠FIR滤波器， 图示滤波器的单位脉冲响应和 损耗函数。 要求通带最大衰减为0.1 dB， 阻带最小衰减为30 dB， 过渡带宽度为0.02 rad。 画出实现抽取器的多相结构， 并求出多相实现时各子滤波器的单位脉冲响应。,题4解图,解： 由式（7.2.2）知道，
16、 抗混叠FIR滤波器的阻带截止频率为s=/D=/5 rad， 要求过渡带宽度为 0.02 rad， 所以， 通带截止频率为/50.02=0.18 rad。 本题已经给出： 通带最大衰减为0.1 dB， 阻带最小衰减为30 dB。 调用MATLAB滤波器设计函数remezord和remez设计该滤波器的程序为ex805.m。 运行程序得到的滤波器单位脉冲响应h(n)和损耗函数如题5解图所示。 请读者仿照教材中图8.6.10画出多相实现结构。（说明： 调用remezord估算得到滤波器长N=M+1=182， 取5的整数倍N=537， 程序中阶数M=M+3。）,本题所给滤波器过渡带太窄， 使h(n)
17、的长度为N=185， 其列表篇幅太大， 所以仅给出h(n)的波形图。 请读者运行程序直接在MATLAB命令窗中查看h(n)的数据， 并根据下式确定多相滤波器实现结构中的5个多相滤波器系数:pk(n)=h(k+nD) k=0, 1, 2, 3, 4; n=0, 1, 2, , 36程序ex805.m清单如下：,% 程序ex805.m % 调用remez函数设计按因子D=5的抽取器的抗混叠滤波器 f=0.18， 1/5； %对归一化边界频率 m=1， 0； rp=0.1； rs=30； dat1=(10(rp/20)1)/(10(rp/20)+1)；dat2=10(rs/20); rip=dat1
18、， dat2； M， fo， mo， w=remezord(f， m， rip)； M=M+3； hn=remez(M， fo， mo， w)； %以下为绘图部分（省略） 运行程序， 绘图如题5*解图所示。,题5*解图,6*. 设计一个内插器， 要求内插因子I=2。 用remez函数设计镜像FIR滤波器， 图示滤波器的单位脉冲响应和损耗函数。 要求通带最大衰减为0.1 dB， 阻带最小衰减为30 dB，过渡带宽度为0.05 rad。 画出实现内插器的多相结构， 并求出多相实现时各子滤波器的单位脉冲响应。,解： 已知内插因子I=2， 所以镜像FIR滤波器的阻带截止频率为s=/I=/2 rad，
19、要求过渡带宽度为0.05 rad， 所以， 通带截止频率为/20.05=0.45 rad。 本题已经给出： 通带最大衰减为0.1 dB，阻带最小衰减为30 dB。 调用MATLAB信号处理工具箱函数remezord和remez设计该滤波器的程序为ex806.m。 运行程序得到的滤波器单位脉冲响应h(n)和损耗函数如题6解图所示。 请读者仿照教材中图8.6.9画出多相实现结构。（说明： 调用remezord估算得到滤波器长度N=M+1=73， 取I的整数倍N=237=74， 程序中阶数用M=M+1调整。）,滤波器单位脉冲响应h(n)数据为 h(0)= 1.250515e004 =h(73) h(
20、1)= 8.034005e003 =h(72) h(2)= 8.680925e003 =h(71) h(3)= 5.520558e004 =h(70) h(4)= 4.859884e003 =h(69) h(5)= 8.774395e004 =h(68) h(6)= 4.896185e003 =h(67) h(7)= 1.961462e003 =h(66) h(8)= 5.290437e003 =h(65) h(9)=3.140232e003 =h(64),h(10)= 5.717291e003 =h(63)h(11)= 4.534276e003 =h(62)h(12)= 6.057635e0
21、03 =h(61)h(13)= 6.217276e003 =h(60)h(14)= 6.256470e003 =h(59)h(15)= 8.190073e003 =h(58)h(16)= 6.193316e003 =h(57)h(17)= 1.053106e002 =h(56)h(18)= 5.803015e003 =h(55)h(19)= 1.328382e002 =h(54)h(20)= 4.965732e003 =h(53),h(21)= 1.651721e002 =h(52)h(22)= 3.536904e003 =h(51)h(23)= 2.036040e002 =h(50)h(24
22、)=1.287036e003 =h(49)h(25)=2.503651e002 =h(48)h(26)=2.139705e003 =h(47)h(27)= 3.097361e002 =h(46)h(28)= 7.432755e003 =h(45)h(29)=3.911230e002 =h(44)h(30)=1.608449e002 =h(43)h(31)= 5.180676e002 =h(42),h(32)= 3.223045e002 =h(41)h(33)=7.689092e002 =h(40)h(34)=7.324125e002 =h(39)h(35)= 1.642101e001 =h(3
23、8)h(36)= 4.345863e001 =h(37)多相滤波器实现结构中的两个多相滤波器系数为p0(n)=h(nI)=h(2n) n=0, 1, 2, , 36p1(n)=h(1+nI)=h(2n+1) n=0, 1, 2, , 36程序ex806.m清单如下： % 程序ex806.m,% 调用remez函数设计按因子I=2的内插器的镜像滤波器f=0.45， 0.5； %对归一化边界频率m=1， 0； rp=0.1； rs=30； dat1=(10(rp/20)1)/(10(rp/20)+1)；dat2=10(rs/20)； rip=dat1， dat2； M， fo， mo， W=rem
24、ezord(f， m， rip)；M=M+1； hn=remez(M， fo， mo， W)； n=0： M； y=n； hn； Mn； fprintf(h(%d)= %9e =h(%d)n， y) %按照对称性列表显示h(n)%以下为绘图检验和h(n)数据显示部分（省略）,题6*解图,运行程序， 绘图如题6*解图所示。,7*. 设计一个按因子2/5降低采样率的采样率转换器， 画出系统原理方框图。 要求其中的FIR低通滤波器过渡带 宽为0.04 rad， 通带最大衰减为1 dB， 阻带最小衰减为 30 dB。 设计FIR低通滤波器的单位脉冲响应， 并画出一种高效实现结构。 解： 按因子2/5降
25、低采样率的采样率转换器的原理方框图如题7解图(一)所示， 图中I=2， D=5。,题7解图(一),根据题目要求可知， 采样率转换器中FIR低通滤波器的技术指标应为： 阻带截止频率s=min/I, /D=/5 rad； 通带截止频率p=s0.04 =0.16 rad； 通带最大衰减p=1 dB， 阻带最小衰减s=30 dB。,滤波器设计程序ex807.m如下： % 程序ex807.m % 调用remez函数设计按因子2/5的采样率转换器中FIR滤波器 f=0.16， 0.2； %对归一化边界频率 m=1， 0； rp=1； rs=30； dat1=(10(rp/20)-1)/(10(rp/20)
26、+1);dat2=10(rs/20); rip=dat1， dat2； M，fo，mo，W=remezord(f， m， rip)；%M=M+2； hn=remez(M， fo， mo， W)； n=0： M； y=n； hn； M-n； fprintf(h(%d)= %9e =h(%d)n， y) %以下为绘图检验部分(省略),运行程序ex807.m， 得到滤波器单位脉冲响应h(n)数据如下： h(0)=1.012427e002 =h(57)h(1)= 1.438004e002 =h(56)h(2)= 1.281009e002 =h(55)h(3)= 1.256566e002 =h(54)h
27、(4)= 1.112633e002 =h(53)h(5)= 7.474111e003 =h(52)h(6)= 1.833721e003 =h(51)h(7)=4.591903e003 =h(50)h(8)=1.004506e002 =h(49),h(9)=1.272172e 002 =h(48) h(10)= 1.139736e 002 =h(47) h(11)= 5.945673e 003 =h(46) h(12)= 2.482502e 003 =h(45) h(13)= 1.155711e 002 =h(44) h(14)= 1.836975e 002 =h(43) h(15)= 2.02
28、0830e 002 =h(42) h(16)= 1.549946e 002 =h(41) h(17)= 4.443471e 003 =h(40) h(18)= 1.063976e 002 =h(39) h(19)= 2.570069e 002 =h(38),h(20)=3.570212e002 =h(37)h(21)=3.593645e002 =h(36)h(22)=2.321793e002 =h(35)h(23)= 2.989859e003 =h(34)h(24)= 4.026557e002 =h(33)h(25)= 8.328551e002 =h(32)h(26)= 1.250531e00
29、1 =h(31)h(27)= 1.581392e001 =h(30)h(28)= 1.763757e001 =h(29),因为DI， 所以， 只要将题7解图（一）中的滤波器h(n)和抽取器用教材第256页中图8.6.2替换， 就得到直接型FIR滤波器高效结构。 由于h(n)较长， 结构图占用篇幅太大， 所以具体结构图请读者练习画出。 当然也可以用多相结构实现， 这时要求滤波器总长度N是抽取因子5的整数倍。 但程序ex807.m运行结果N=M+1=58， 所以， 应当保留程序中第7行后面的阶数修改语句“M=M+2”， 运行得到新的h(n)， 其长度N=60=512， 画出多相结构如题7解图(二)
30、所示。,题7解图（二）,图中5个多相滤波器系数按照如下公式确定： pk(n)=h(k+5n) k=0, 1, 2, 3, 4; n=0, 1, 2, , 11 请读者运行程序并确定pk(n)的具体数据。 8*. 假设信号x(n)是以奈奎斯特采样频率对模拟信号xa(t)的采样序列， 采样频率Fx =10 kHz。 现在为了减少数据量， 只保留0f3 kHz的低频信息， 希望尽可能降低采样频率， 请设计采样率转换器。 要求经过采样率转换器后， 在频带0f2.8 kHz中频谱失真不大于1 dB， 频谱混叠不超过1。,（1） 确定满足要求的最低采样频率Fy和相应的采样率转换因子;（2） 画出采样率转换
31、器原理方框图;（3） 确定采样率转换器中FIR低通滤波器的技术指标， 用等波纹最佳逼近法设计FIR低通滤波器， 画出滤波器的单位脉冲响应及其损耗函数曲线， 并标出指标参数（通带截止频率、 阻带截止频率、 通带最大衰减和阻带最小衰减）；（4） 求出多相实现结构中子滤波器的单位脉冲响应， 并列表显示或打印。,解： （1） 根据时域采样定理， 满足要求的最低采样频率Fy=23 kHz=6 kHz， 采样率转换因子为Fy/Fx=6/10=3/5。（2） 采样率转换器原理方框图如题8解图(一)所示。 图中， I=3， D=5。,题8解图(一),（3） 根据题目要求可知， 采样率转换器中FIR低通滤波器的
32、技术指标应为： 通带截止频率p=22800/(IFx)= 22800/30000=28/150 rad（请读者注意， 滤波器工作 在内插后及抽去之前， 所以应当用采样频率IFx计算2800 Hz对应的数字频率）， 阻带截止频率s=min/I, /D= /5 rad； 根据要求“在频带0f2.9 kHz中频谱失真不大于1 dB”得到： 通带最大衰减p=1 dB， 根据要求“频谱混叠不超过1”可知， 阻带最小衰减s=20 lg(1%)=40 dB。 FIR低通滤波器的设计程序为ex808.m， 运行程序画出滤波器的单位脉冲响应及其损耗函数曲线如题8解图(二)所示。,题8解图(二),说明： 调用re
33、mezord估算得到滤波器长度N=M+1=214， 为了满足多相实现条件， 取D的整数倍N=543=215， 程序中阶数用M=M+1调整。 （4） 按照下式构造多相实现结构中子滤波器的单位脉冲响应： pk(n)=h(k+5n) k=0, 1, 2, 3, 4; n=0, 1, 2, , 42 由于本题给定的指标太高， 使滤波器阶数高达214， 列出多相实现结构中子滤波器的单位脉冲响应序列的数据篇幅太大， 所以， 请读者运行程序显示h(n)数据， 并确定5个子滤波器的单位脉冲响应序列pk(n)的具体数据。,本题求解程序ex808.m清单如下： % 程序ex808.m % 调用remez函数设计按因子3/5的采样率转换器中FIR滤波器 f=28/150， 1/5； %对归一化边界频率 m=1， 0； rp=1； rs=40； dat1=(10(rp/20)1)/(10(rp/20)+1)； dat2=10(rs/20)； rip=dat1， dat2；,M， fo， mo， W=remezord(f， m， rip)；M=M+1； %阶数M+1， 使N=M+1满足5的整数倍要求hn=remez(M， fo， mo， W)； n=0： M； y=n； hn； Mn； fprintf(h(%d)= %9e =h(%d)n， y) %按照对称性列表显示h(n)%以下为绘图检验部分（省略）,
