第8章 多采样率数字信号处理 第第8章章 多采样率数字信号处理多采样率数字信号处理8.1 引言引言 8.2 信号的整数倍抽取信号的整数倍抽取 8.3 信号的整数倍内插信号的整数倍内插 8.4 按有理数因子按有理数因子I/D的采样率转换的采样率转换 8.5 整数倍抽取和内插在数字语音系统中的应用整数倍抽取和内插在数字语音系统中的应用 8.6 采样率转换滤波器的高效实现方法采样率转换滤波器的高效实现方法 8.7 采样率转换器的采样率转换器的MATLAB实现实现 习题与上机题习题与上机题第8章 多采样率数字信号处理 8.1 引引 言言前面所讨论的信号处理的各种方法都是把采样率Fs视为固定值， 即在一个数字系统中只有一个采样频率。 但在实际系统中， 经常会遇到采样率的转换问题， 即要求一个数字系统能工作在“多采样率”状态。 例如： 第8章 多采样率数字信号处理 （1） 在数字电视系统中， 图像采集系统一般按444标准或422标准采集数字电视信号， 再根据不同的电视质量要求， 将其转换成其它标准的数字信号（如422， 411， 211等标准）进行处理、 传输。 这就要求数字电视演播室系统工作在多采样率状态。 （422标准的含义是“亮度信号Y的采样率: 红色差信号R-Y的采样率:蓝色差信号B-Y的采样率=422”， 其他标准以此类推。）第8章 多采样率数字信号处理 （2） 在数字电话系统中， 传输的信号既有语音信号， 又有传真信号， 甚至有视频信号， 这些信号的带宽相差甚远。 所以， 该系统应具有多采样率功能， 并根据所传输的信号自动完成采样率转换。 （3） 对一个非平稳随机信号（如语音信号）作谱分析或编码时， 对不同的信号段， 可根据其频率成分的不同而采用不同的采样率， 以达到既满足采样定理， 又最大限度地减少数据量的目的。 （4） 如果以高采样率采集的数据存在冗余， 这时就希望在该数字信号的基础上降低采样速率， 剔除冗余， 减少数据量， 以便存储、 处理与传输。第8章 多采样率数字信号处理 以上所列举的几个方面都是希望能对采样率进行转换， 或要求数字系统工作在多采样率状态。 近年来， 建立在采样率转换基础上的“多采样率数字信号处理”已成为数字信号处理学科的主要内容之一。 一般认为， 在满足采样定理的前提下， 首先将以采样率F1采集的数字信号进行D/A转换, 变成模拟信号， 再按采样率F2进行A/D变换， 从而实现从F1到F2的采样率转换。 但这样较麻烦， 且易使信号受到损伤， 所以实际上改变采样率是在数字域实现的。 根据采样率转换理论， 对采样后的数字信号x(n)直接进行采样率转换， 以得到最新采样率下的采样数据。 第8章 多采样率数字信号处理 采样率转换通常分为“抽取（Decimation）”和“插值（Interpolation）”。 抽取是降低采样率以去掉多余数据的过程， 而插值则是提高采样率以增加数据的过程。 本章先讨论抽取和插值的一般概念， 然后讨论其几种基本的实现方法。 本章所涉及的内容也是语音及图像数据压缩新技术子带编码的重要理论基础。 第8章 多采样率数字信号处理 8.2 信号的整数倍抽取信号的整数倍抽取设x(n1T1)是连续信号xa(t)的采样序列， 采样率F1=1/T1(Hz)， T1称为采样间隔， 单位为秒， 即如果希望将采样率降低到原来的1/D， D为大于1的整数， 称为抽取因子。 最简单的方法是对x(n1T1)每D点抽取1点， 抽取的样点依次组成新序列y(n2T2)。 y(n2T2)的采样间隔为T2， 采样率为F2 = 1/T2(Hz)， T2与T1的关系为（）第8章 多采样率数字信号处理 为了后面叙述方便， 我们将上述的抽取系统用图8.2.1(a)表示， 图中符号表示采样率降低为原来的1/D（D为Decimation的第一个字母， 表示抽取）。 x(n1T1)和y(n2T2)分别如图8.2.1（b）和（c）所示。 图中n1和n2分别表示x(n1T1)和x(n2T2)序列的序号， 于是有（8.2.3）当n1=n2D时， y(n2T2)=x(n1T1)。 （8.2.2）D 第8章 多采样率数字信号处理 图8.2.1 数字信号的时域抽取示意图 第8章 多采样率数字信号处理 上面在时域讨论了整数倍抽取的概念。 抽取看起来好像很简单， 只要每隔D1个抽取一个就可以了， 但抽取降低了采样频率， 会引起频谱混叠现象。 下面讨论抽取过程中可能出现的频谱混叠及改进措施。 如果x(n1T1)是连续信号xa(t)的采样信号， 则xa(t)和x(n1T1)的傅里叶变换Xa(j)和将分别是 （） （） 第8章 多采样率数字信号处理 其中， =2f rad/s， f为模拟频率变量; 1为数字频率,（）由（）式有 （） 式中， sa1=2/T1 rad/s, 亦称为采样频率。 第8章 多采样率数字信号处理 为了对抽样前后的频谱进行比较， 作图时均以模拟角频率为自变量（横坐标）， 为此按（）式将写成的函数为 （） 因为这里xa(t)是一般的非周期连续函数， 所以Xa(j)也是模拟频率的非周期函数， 如图（a）所示。 第8章 多采样率数字信号处理 而x(n1T1)的傅里叶变换为连续频率1的周期函数。 在满足采样定理时，的频谱在sa1/2, sq1/2上与Xa(j)相似（差一个比例常数1/T1）， 且无混叠现象， 如图8.2.2（b）所示。 但如果将采样率降低到原来的1/D， 即T2=DT1， 当D=4时， 得到y(n2T2)及其频谱如图8.2.3所示（实际上， 应为图中各重复谱的叠加曲线）。 图中， y(n2T2)为对x(n1T1)抽取的结果， 为y(n2T2)的傅里叶变换。 的周期sa2=2/T2=2/DT1=sa1/D。 这就是说， 的周期是周期的1/D。 第8章 多采样率数字信号处理 图8.2.2 xa(t)与x(n1T1)及其频谱图第8章 多采样率数字信号处理 图8.2.3 抽取引起的频谱混叠现象第8章 多采样率数字信号处理 由图8.2.3可见， 是有混叠的， 无法从y(n2T2)中恢复出x(n1T1)来。 所以随意对x(n1T1)进行抽取是不行的。 只有在抽取后仍能满足采样定理时才能恢复出原来的信号xa(t)， 否则就必须另外采取措施。 通常采取的措施是抗混叠滤波。 所谓抗混叠滤波， 就是在抽取之前先对信号进行低通滤波， 把信号的频带限制在sa2/2以下。 这种抽取系统框图如图8.2.4所示。 图中h(n1T1)为抗混叠滤波器， 它的输出v(n1T1)的最高频率已被h(n1T1)限制在sa2/2=sa1/(2D)以下。 即抗混叠滤波器的阻带截止频率为sa1/(2D)， 对应的数字阻带截止频率为第8章 多采样率数字信号处理 图8.2.4 带有抗混叠滤波器的抽取系统框图第8章 多采样率数字信号处理 所以， 在理想情况下， 抗混叠低通滤波器h(n1T1)的频率响应H(ej)由下式给出: （8.2.9） 第8章 多采样率数字信号处理 图8.2.4中各点的信号在时域和频域中的示意图如图8.2.5所示。 这种办法虽然把x(n1T1)中的高频部分损失掉了， 但由于抽取后避免了混叠， 所以在中完好无损地保留了中的低频部分， 可以从中恢复出的低频部分。 第8章 多采样率数字信号处理 图8.2.5 抽取前后信号的时域和频域示意图第8章 多采样率数字信号处理 为了进一步搞清x(n1T1)经 前后的频谱关系， 对信号的抽取过程进行等效的数学描述， 如图8.2.6所示， 以便于进行频谱分析， 找出与之间的关系。 在抽取前先令x(n1T1)乘以周期序列， 即（8.2.10） 其中， 定义如下： Ddef第8章 多采样率数字信号处理 D然后对 进行， 得到y(n2T2)。 的离散傅里叶级数（DFS）系数为（8.2.11） 于是的DFS展开式为（8.2.12） 第8章 多采样率数字信号处理 将（8.2.12）式代入（8.2.10）式, 得（8.2.13）对进行抽取， 每隔D-1个点抽取一个样值， 所取得抽样点均在的点上。 把这样抽取的结果作为y(n2T2)， 如图8.2.6所示（D=4）。 第8章 多采样率数字信号处理 图8.2.6 抽取过程的等效数学描述与直接抽取波形第8章 多采样率数字信号处理 下面推导与的关系： 当n2D=n1时， ， 而当n2Dn1时， =0， 所以 （8.2.14）第8章 多采样率数字信号处理 令 , 则（8.2.15）式中， 。 所以有（省去z的下标）（8.2.16） （8.2.14）式就是与的关系， 即是的D个平移样本之和， 相邻的平移样本在频率轴1上相差2/D， 在模拟频率轴上相差2/(DT1)=sa1/D=sa2， 如图8.2.7和图8.2.8所示。第8章 多采样率数字信号处理 为了更清楚地说明上述关系， 将由求得y(n2T2)频谱的过程用图和图表示出来。 图给出了抽取后产生混叠的情况， 即中的最高频率csa1/(2D)=sa2/2。 这里sa1为的周期， sa2为的周期。 图给出了抽取后不产生混叠的情况， 即中的最高频率csa2/2时， 抽取前后的频谱关系示意图第8章 多采样率数字信号处理 图8.2.8 c2c， 所以允许有一定的过渡带， 可用线性相位FIR滤波器实现。 根据其功能， 将h(n2T2)称为镜像滤波器。 第8章 多采样率数字信号处理 图8.3.6 零值内插前后的时域信号及其频谱 第8章 多采样率数字信号处理 图8.3.7 镜像滤波器的理想幅频特性第8章 多采样率数字信号处理 将理想镜像滤波器的阻带截止频率换算成数字频率为所以， 理想情况下， 镜像滤波器h(n2T2)的频率响应特性为 （8.3.4）第8章 多采样率数字信号处理 式中， C为定标系数。 因此输出频谱为 （8.3.5）定标系数C的作用是， 在m=0, I, 2I, 3I, 时， 确保输出序列y(m)=x(m/I)。 为了计算简单， 取m=0来求解C的值。 第8章 多采样率数字信号处理 因为2=1/I， 所以由此得出， 定标系数C=I。 第8章 多采样率数字信号处理 3 内插器的输入、内插器的输入、 输出关系输出关系1） 时域输入、 输出关系由图有因为（）第8章 多采样率数字信号处理 所以上式就是内插器时域输入、 输出关系。 （）第8章 多采样率数字信号处理 由（8.3.3）式知道， 所以在复频域分析图8.3.2时， 其输入x(n1T1)的Z变换X(z1)与输出y(n2T2)的Z变换Y(z2)的关系推导如下： （8.3.10）（8.3.9）（）2） 频域输入、 输出关系第8章 多采样率数字信号处理 其中（8.3.11） （8.3.12） (8.3.12)式中所有变量都为z2， 所以可去掉下标， 得到：(8.3.13)为的整数倍即时第8章 多采样率数字信号处理 8.4 按有理数因子按有理数因子I/D的采样率转换的采样率转换在按整数因子I内插和整数因子D抽取的基础上， 本节介绍按有理数因子I/D采样率转换的一般原理。 显然， 可以用图所示方案实现有理数因子I/D采样率转换。 第8章 多采样率数字信号处理 图8.4.1 按有理数因子I/D的采样率转换方法第8章 多采样率数字信号处理 首先对输入序列x(n)按整数因子I内插， 然后再对内插器的输出序列按整数因子D抽取， 达到按有理数因子I/D的采样率转换。 应当注意， 先内插后抽取, 才能最大限度地保留输入序列的频谱成分。 用Fx=1/Tx和FY=1/Ty分别表示输入序列x(n)和输出序列y(m)的采样频率， 则Fy=(I/D)Fx。 另外， 图中镜像滤波器hI(l)和抗混叠滤波器hD(l)级联， 而且工作在相同的采样频率IFx， 因此完全可以将它们合成为一个等效滤波器h(l)， 得到按有理数因子I/D采样率转换的实用原理方框图， 如图所示。 第8章 多采样率数字信号处理 图8.4.2 按有理数因子I/D采样率转换的实用原理方框图第8章 多采样率数字信号处理 如前所述， 理想情况下， hI(l)和hD(l)均为理想低通滤波器， 所以等效滤波器h(l)仍是理想低通滤波器， 其等效带宽应当是hI(l)和hD(l)中最小的带宽。 h(l)的频率响应为（8.4.1） 第8章 多采样率数字信号处理 现在推导图中输出序列y(m)的时域表达式。 零值内插器的输出序列为（）线性滤波器输出序列为（）其它第8章 多采样率数字信号处理 整数因子D抽取器最后输出序列为y(m)， 其时域表达式为 （8.4.4）如果线性滤波器用FIR滤波器实现， 则可以根据式（8.4.4）计算输出序列y(m)。 除了前面介绍的采样率变换技术， 在实际工作中还会遇到任意因子采样率转换（Fy/Fx为任意有限数或可能随机变化）。 有兴趣的读者请参考文献30。 第8章 多采样率数字信号处理 为了下面叙述方便， 首先说明本节对信号时域和频域的表示方法和描述符号。 设x(t)为模拟信号， x(nT1)表示对x(t)的采样序列， y(mT2)是对x(nT1)进行采样率转换（内插或抽取）后的序列。 并定义 8.5 整数倍抽取和内插在数字语音系统中的应用整数倍抽取和内插在数字语音系统中的应用第8章 多采样率数字信号处理 其中， 数字频率与模拟频率的关系为1=1, 2=2。 x(nT1)的采样频率记为Fsa1=1/T1 Hz， y(nT2)的采样频率记为Fsa2=1/T2 Hz， 相应的采样角频率记为sa1=2Fsa1=2T1 rad/s， sa2=2Fsa2=2/T2 rad/s。 为了通过观察比较x(t)、 x(nT1)和y(mT2)的频谱关系， 理解采样率转换在数字语音系统中的应用原理， 本节全部以模拟角频率为自变量（横坐标）， 并采用上面定义的符号， 来绘制应用系统中各信号的频谱曲线。 第8章 多采样率数字信号处理 8.5.1 数字语音系统中的信号采样过程及数字语音系统中的信号采样过程及其存在的问题其存在的问题在数字语音系统中， 语音信号的采样过程如图所示。 图中， x(t)为模拟信号， 其有用频谱分布范围为fh, fh， fh表示x(t)中有用频率成分的最高频率。 信号中一般含有干扰噪声， 其频带宽度远大于fh。 x(t)及其幅频特性|X(j)|如图（b）所示。 下面以电话系统中的数字语音系统为例， 讨论图（a）所示的基本采集系统中存在的技术问题。第8章 多采样率数字信号处理 图8.5.1 语音信号的一般采样过程示意图第8章 多采样率数字信号处理 在电话系统中， 一般要保证4 kHz的音频带宽， 即取fh=4 kHz。 但送话器发出的信号x(t)的带宽比fh大很多。 因此， 在A/D变换之前要对其进行模拟预滤波， 以防止采样后发生频谱混叠失真。 为了使信号采集数据量尽量小， 取采样频率Fs=2fh=8 kHz。 这时要求低通模拟滤波器h(t)的幅频响应特性|H(j)|如图（c）所示。 预滤波后的信号v(t)及其采样序列v(nT)和相应的频谱分别如图（d）、 （e）所示。 第8章 多采样率数字信号处理 上述基本集系统对x(t)进行A/D变换的困难在于对预滤波器h(t)的技术要求太高（要求过渡带宽度为0， 用理想低通滤波器）， 因而是难以设计与实现的。 显然， 在接收端D/A变换过程中同样会遇到此问题。 如果简单地将采样率提高， 如取Fs=16 kHz， 则预滤波器就容易实现（允许有4 kHz的过渡带）， 但使采集信号的数据量加大1倍， 传输带宽也加大1倍。 下面讨论如何采用整数因子抽取与整数因子内插来解决该问题， 而不增加数据量。第8章 多采样率数字信号处理 8.5.2 数字语音系统中改进的数字语音系统中改进的A/D转换方案转换方案为了降低对模拟预滤波器的技术要求， 采用如图8.5.2（a）所示的改进方案。 先用较高的采样率进行采样， 如采样率Fsa1=1/T1=16 kHz， 经过A/D后， 再按因子D=2抽取， 把采样率降至8 kHz。 这时， 模拟预滤波器g(t)的过渡带为4f12 kHz， 如图8.5.2（c）所示。 这样的预滤波器会导致采样信号w(nT1)的频谱在412 kHz的频带中发生混叠， 如图8.5.2（e）所示。 但这部分混叠在抽取前用数字滤波器h(nT1)滤掉了。第8章 多采样率数字信号处理 数字滤波器h(nT1)的幅频特性|如图8.5.2（f）所示。 这样， 模拟预滤波器就容易设计和实现了。 现在把问题转移到设计和实现技术要求很高的数字滤波器h(nT1)上了， 这就是解决问题的关键技术。 数字滤波器可用FIR结构， 容易设计成线性相位和陡峭的通带边缘特性。 这种方案最终并未增加信号数据量。 )e (1j THW第8章 多采样率数字信号处理 图8.5.2 数字语音系统中改进的A/D转换方案及各点信号波形与频谱示意图第8章 多采样率数字信号处理 8.5.3 接收端接收端D/A转换器的改进方案转换器的改进方案设数字信号序列y(mT2)传送到接收端后变成， 若不考虑信道噪声， 则其频谱与图8.5.2（h）相同。 要将恢复为模拟信号， 若采用基本方案， 先将经D/A转换器， 再进行模拟低通滤波， 得到。 这种方案同样会对模拟恢复低通滤波器提出难以实现的技术要求。 为了解决这一难题， 可采用如图8.5.3所示D/A转换器的改进方案。 该方案的思路是， 采用整数因子内插， 将模拟恢复低通滤波器的设计与实现的困难转移到设计滤除镜像频谱的高性能数字滤波器来解决。 具体实现原理如下。 第8章 多采样率数字信号处理 图8.5.3 D/A转换器的改进方案第8章 多采样率数字信号处理 设输入数字信号 如图8.5.4（a）所示（与图8.5.2（h）相同）。 经内插后将采样率提高2倍， 滤波器的输出为， 假定可设计成陡峭通带边缘特性， 则的时域和频域波形如图8.5.4（b）所示。 对进行D/A变换， 得到:(8.5.1)及其幅频特性如图8.5.5所示。 应当说明， 这种D/A转换器难以实现， 实际中常用零阶保持型D/A转换器代替， 但其频响特性不理想， 会引入幅频失真。 这种失真可在数字域进行预处理补偿。 第8章 多采样率数字信号处理 图8.5.4 和时域和频域示意图第8章 多采样率数字信号处理 图8.5.5 时域和频域示意图第8章 多采样率数字信号处理 对进行模拟低通滤波， 这时要求模拟低通滤波器的通带边缘频率为p=/(2T1)， 过渡带为/(2T1)|3/(2T1)， 阻带为3/(2T1)|。 的幅频特性曲线如图8.5.6所示， 当然， 过渡带上的频响曲线可以不是直线。 的输出则为模拟信号。 由于过渡带较宽， 所以模拟低通滤波器的设计与实现较容易。 我们希望恢复的信号就是， 其时域和频域示意图如图8.5.7所示。 第8章 多采样率数字信号处理 图8.5.6 的幅频特性曲线 第8章 多采样率数字信号处理 图8.5.7 恢复的模拟信号 及其频谱示意图 第8章 多采样率数字信号处理 8.6 采样率转换滤波器的高效实现方法采样率转换滤波器的高效实现方法8.6.1 直接型直接型FIR滤波器结构滤波器结构1 整数倍抽取器的整数倍抽取器的FIR直接实现直接实现整数（D）倍抽取器框图如图8.2.4所示。 抗混叠低通滤波器用FIR结构时， 抽取器的时域输入、 输出关系为（设h(n1T1)长度为N）（8.6.1）（8.6.2）第8章 多采样率数字信号处理 如果滤波器用FIR直接型结构， 则该抽取器的实现网络结构如图（a）所示。 经滤波卷积运算得出v(n1T1）， 最后将v(n1T1)每隔D1个取一个作为输出y(n2T2)， 即v(n1T1)中有（D1）/D的样值都被舍弃了。 所以这种结构是一种低效实现结构， 而且要求计算每一个v(n1T1)的N次乘法和N1次加法在一个T1时间内完成。 第8章 多采样率数字信号处理 图8.6.1 按整数因子D抽取系统的直接型FIR滤波器结构第8章 多采样率数字信号处理 为了得到相应的高效FIR直接实现， 对图8.6.1（a）进行等效变换。 显然， 将图8.6.1（a）中的移在N条乘法器支路中的乘法器之前， 如图8.6.1（b）所示， 所得y(n2T2)与原结构输出相同， 即图8.6.1（a）与图8.6.1（b）是等效的。 图8.6.1（b）中各条支路里的同时在n1=n2D时开通， 例如D=4， N=11， n2=5时， n1=20， 第0条支路通过的是x(20T1)， 第1条支路通过的是x(19T1)， 最下面的第N1条支路通过的是x(10T1)。 DD第8章 多采样率数字信号处理 此刻开始计算N个支路的N次乘法和最后的N1次加法， 得到一个输出样值： y(n2T2)=y(5T2)=y(5DT1)=y(20T1)。 由于在x(21T1)到来之前所有的同时关闭， 直到n2=6时， 即n1=n2D=24时， N个才又同时开通， 分别让x(24T1)， x(23T1)， ， x(11T1)通过， 开始计算下一个输出序列样值y(6T2)。 所以， 改进后的实现结构将乘法运算移到低采样率一侧， 使乘法运算速度要求降低到原来的1/D， 即原来要在一个T1时间内完成的运算， 现在只要在DT1时间之内完成就可以了。 当然， 也使计算量减少到原来的1/D。 故称之为高效结构。 DD第8章 多采样率数字信号处理 应当说明， 图8.6.1（b）中将放在h(0)， h(T1)， ， h(N1)T1之前， 减少了运算量， 但这并不是把抗混叠滤波放到了抽取之后， 而是与原来的滤波作用等效。 对此作如下解释： 滤波和抽取的作用次序在FIR实现结构中体现在滤波器输入端及延迟链上所加的信号序列， 如果所加信号是抽取以前的信号， 则是先滤波后抽取， 反之是先抽取后滤波。 图8.6.1（b）中， 所有均安排在延迟链之后， 即滤波器延迟链上各点的信号仍然是原序列x(n1T1)， x(n11)T1， ， x(n1N+1)T1， 而不是抽取后的信号。DD第8章 多采样率数字信号处理 每当开通时， 进入左侧的信号是未抽取的原信号， 即输出的y(n2T2)与图8.6.1（a）中抽选的y(n2T2)相同， 而两次开通之间所阻挡的信号恰好就是图8.6.1（a）中将来要舍弃的部分， 所以计算结果是正确的。 但绝对不能将提前到延迟链之前， 那样才是真正的先抽取后滤波器， 会产生严重的混叠现象。 DD第8章 多采样率数字信号处理 由于常常希望把FIR滤波器设计成线性相位的， 因而实现设计滤波器时可用FIR线性相位结构， 这样又可以使乘法计算量减少一半。 根据线性相位时域特性可画出抽取器FIR结构的线性相位形式如图8.6.2所示。 第8章 多采样率数字信号处理 图8.6.2 抽取器FIR结构的线性相位形式(N=11)第8章 多采样率数字信号处理 2 整数倍内插器的整数倍内插器的FIR直接实现直接实现整数倍内插系统框图如图所示。 镜像滤波器h(n2T2)采用FIR结构时， I倍内插器的FIR直接实现结构如图所示。 图中乘法是在高采样率一侧进行的， 不是高效结构， 应设法将乘法运算移到低采样率一侧以减少计算量。 可采用以下方法进行网络等效变换， 得出相应的高效结构。 先将FIR滤波网络部分进行转置， 得到图所示的FIR转置型结构， 再用其代替图中的FIR滤波网络， 得到图所示的内插系统直接实现。 第8章 多采样率数字信号处理 图8.6.3 按整数因子I内插系统的直接型FIR滤波器结构 第8章 多采样率数字信号处理 图8.6.4 FIR转置型结构 第8章 多采样率数字信号处理 由整数倍抽取系统的FIR直接实现的等效变换概念可知， 图8.6.5中先零值内插后分支相乘与先分支相乘后零值内插等效。 因此， 可将图8.6.5中的分别移到FIR网络的各支路的乘法器之后， 可得到图8.6.6所示的内插系统直接实现高效结构。 由于延时链上所加的仍然是内插后的信号， 因而等效变换后的高效结构仍是先内插后滤波。 第8章 多采样率数字信号处理 图8.6.5 滤波网络转置后的内插器的第8章 多采样率数字信号处理 图8.6.6 内插器FIR直接实现的高效结构 第8章 多采样率数字信号处理 当满足线性相位条件h(n2T2)=h(N1n2)T2)时， 可用线性相位结构实现， 将乘法次数再减少一半。 取N=9， 画出内插器的线性相位FIR直接高效实现, 如图所示。 第8章 多采样率数字信号处理 图8.6.7 内插器的线性相位FIR直接实现结构第8章 多采样率数字信号处理 观察图和图（b）可发现一个有趣的规律: 图所示的按整数因子I内插系统的高效FIR滤波器结构与图（b）所示的按整数因子D抽取系统的高效FIR滤波器结构互为转置关系。 这种关系有助于对简化整数因子抽取系统和整数因子内插系统的高效FIR滤波器结构的研究， 在多相实现结构的讨论中将用到该关系。 第8章 多采样率数字信号处理 3 按有理数因子按有理数因子I/D的采样率转换系统的高效的采样率转换系统的高效FIR滤滤波器结构波器结构为了叙述方便， 先由图画出按有理数因子I/D采样率转换系统的直接型FIR结构， 如图所示。 第8章 多采样率数字信号处理 图8.6.8 按有理数因子I/D采样率转换系统的直接型FIR滤波器结构第8章 多采样率数字信号处理 按有理数因子I/D采样率转换系统的高效结构一般基于内插系统的高效FIR滤波器结构与抽取系统的高效FIR滤波器结构进行设计。 其指导思想是， 使FIR滤波器运行于最低采样速率。 为此， 当ID时， Fy Fx， 将图8.6.8中的直接型FIR结构与前面的用图8.6.6所示的整数因子I内插器的高效FIR滤波器结构代替即可。 当ID时， FyFx， 将图8.6.8中的直接型FIR结构与后面的 用图8.6.1（b）所示的整数因子D抽取器的高效FIR滤波器结构代替即可。 如果采用线性相位FIR滤波器， 则应当用相应的线性相位FIR滤波器的高效内插结构或高效抽取结构来实现。 ID第8章 多采样率数字信号处理 （） 8.6.2 多相滤波器结构多相滤波器结构可以证明， 图8.6.6所示的按整数因子I内插系统的高效FIR滤波器结构可以用一组较短的多相滤波器组实现。 如果FIR滤波器总长度为M=NI， 则多相滤波器组由I个长度为N=M/I的短滤波器构成， 且I个短滤波器轮流分时工作。 为了证明上述结论， 观察图8.6.3给出的整数因子I内插系统的直接型FIR滤波器结构。 为了下面描述简单， 用x(n)表示x(n1T1)， 用v(m)表示v(n2T2)， y(m)表示y(n2T2)。 输出序列y(m)为第8章 多采样率数字信号处理 零值内插器的输出序列v(m)是在输入序列x(n)的两个相邻样值之间插入I1个零样值得到的， 因此v(m)进入FIR滤波器的M个样值中只有N=M/I个非零值。 所以在任意m时刻， 计算y(m)=h(m)*v(m)时只有N个非零值与h(m)中的N个系数相乘。 由（8.3.2）式知道所以， m=jI时刻， 第8章 多采样率数字信号处理 (8.6.4) m=jI+1时刻， （8.6.4）中v(jIn)右移1位， N个x(n)的非零值与h(n)的对应关系又右移1位， 所以(8.6.5) 第8章 多采样率数字信号处理 当m=jI+I时刻， N个x(n)的值与h(n)的对应关系又重复（8.6.4）式， 只是x(n)又移进1位， 所以(8.6.6) 综上所述， 当m=jI+k, k=0, 1, 2, , I1, j=0, 1, 2, 时, 有（8.6.7） 第8章 多采样率数字信号处理 把式（8.6.7）中的h(k+nI)看做长度N=M/I的子滤波器的单位脉冲响应， 并用pk(n)表示： 这样， 从m=0开始， 整数因子I内插系统的输出序列y(m)计算如下： （8.6.9）（8.6.8）第8章 多采样率数字信号处理 式中， m=jI+k; k=0, 1, 2, , I-1; j=0, 1, 2, 。 显然， 当m=jI+k从0开始增大时， k从0开始以I为周期循环取值； j表示循环周期数。 所以， 实现（）式的多相滤波器结构如图所示 。 I个子滤波器均运行于低采样率Fx下， 且系数少， 计算量小， 所以多相滤波器结构是一种高效结构。 输入端的x(n)每移入一个样值， I个子滤波器分别计算出y(m)的I个样值， 选择电子开关以高采样率Fy=IFx， 依次逆时针循环选取I个子滤波器的输出， 形成输出序列y(m)。 实现了整数因子I内插功能。 第8章 多采样率数字信号处理 图8.6.9 整数因子I内插系统的多相滤波器结构第8章 多采样率数字信号处理 从I个子滤波器的频响特性可以解释“多相滤波器”的含义。 对低通滤波器h(n)按整数因子I抽取得到子滤波器pk(n)。 h(n)是截止频率为/I的理想低通滤波器， 则pk(n)的截止频率必然是， 即I个子滤波器都是全通滤波器， 幅度特性相同， 它们的唯一区别是相位特性不同， 故称为“多相滤波器”结构。 正是由h(n)的I个不同的起始点抽取得到I个子滤波器， 形成了这种多相特性。 根据整数因子I内插器的实现结构与整数因子D抽取器的实现结构互为转置关系的规律， 将图8.6.9 给出的整数因子I内插系统的多相滤波器结构进行转置， 则得到图 8.6.10 所示的整数因子D抽取系统的多相滤波器结构。 第8章 多采样率数字信号处理 图8.6.10 整数因子D抽取系统的多相滤波器结构第8章 多采样率数字信号处理 定义多相滤波器的单位脉冲响应为 式中， N为pk(n)的长度。 一般选择原抗混叠FIR滤波器总长度M=DN， N=M/D。 电子开关以速率Fx逆时针旋转， 从子滤波器p0(n)在m=0时刻开始， 并输出y(0)； 然后电子开关以速率Fx逆时针每旋转一周， 即每次转到子滤波器p0(n)时， 输出端就以速率Fy=Fx/D送出一个y(m)样值。 （8.6.10）第8章 多采样率数字信号处理 下面以N=D=2， M=DN=4为例， 验证图8.6.9所示的抽取系统多相结构的正确性。 首先根据图8.6.1（a）计算出抽取器的正确输出y(m): 假设x(n)为因果信号， 则 （8.6.11） 第8章 多采样率数字信号处理 现在根据图计算多相实现结构的输出y(m)。 图中， 多相子滤波器p0(n)=h(0), h(2), p1(n)=h(1), h(3), 开始时k=0, n=0， 只有x(0)进入p0(n)， p1(n)中无信号， 所以总输出y(0)=p0(0)x(0)=h(0)x(0)。 逆时针旋转开始下一周期: k=D1=1时， 电子开关转到p1(n), x(1)进入p1(n)， p1(n)的输出为p1(0)x(1)=h(1)x(1)； k=0时， 电子开关又转到p0(n)， 此时， x(2)进入p0(n)第一节， 上一周期中进入p0(n)的x(0)移位到p0(n)的第二节， 所以p0(n)的输出为p0(0)x(2)+p0(1)x(0)=h(0)x(2)+h(2)x(0)总的输出y(1)为p0(n)与p1(n)输出之和， 即y(1)=h(0)x(2)+h(2)x(0)+h(1)x(1)第8章 多采样率数字信号处理 同样道理， 可求出下一旋转周期得到的输出为所求y(0)、 y(1)和y(2)与（8.6.11）式相同， 所以， 图8.6.10所给结构是正确的。 第8章 多采样率数字信号处理 需要说明， 在实际采样率转换系统中， 常常会遇到抽取因子和内插因子很大的情况。 例如， 按有理数因子I/D=150/61的采样率转换系统， 从理论上讲， 可以采用多相滤波器结构准确地实现这种采样率转换， 但是实现结构中将需要150个多相滤波器， 而且其工作效率很低。 “多级实现结构”可以很好地解决该问题。 而且， 多级实现结构可以使滤波器总长度大大降低。 多级实现内容请参考文献12， 19。 按有理数因子I/D采样率转换系统还可以采用线性时变滤波器结构实现12。 第8章 多采样率数字信号处理 【例例】 设计一个按因子I=5的内插器， 要求镜像滤波器通带最大衰减为0.1 dB， 阻带最小衰减为30 dB， 过渡带宽度不大于/20。 设计FIR滤波器系数h(n)， 并求出多相滤波器实现结构中的5个多相滤波器系数。 解解 由（）式知道FIR滤波器h(n)的阻带截止频率为/5， 根据题意可知滤波器其他指标参数: 通带截止频率为/5/20=3/20， 通带最大衰减为0.1 dB， 阻带最小衰减为30 dB。 调用remezord函数求得h(n)长度M=47， 为了满足5的整数倍， 取M=50。 调用remez函数求得h(n)如下:第8章 多采样率数字信号处理 h002=h(49) h(1)=002=h(48)h(2)=002=h(47)h(3)=002=h(46)h(4)=002=h(45)h(5)=002=h(44)h(6)=002=h(43)h(7)=002=h(42)h002=h(41)第8章 多采样率数字信号处理 h002=h(40)h002=h(39)h002=h(38)h002=h(37)h(13)=003=h(36)h(14)=002=h(35)h(15)=001=h(34)h(16)=001=h(33)第8章 多采样率数字信号处理 h(17)=001=h(32)h(18)=002=h(31)h002=h(30)h001=h(29)h001=h(28)h001=h(27)h001=h(26)h001=h(25)第8章 多采样率数字信号处理 根据（）式确定多相滤波器实现结构中的5个多相滤波器系数如下:p0(n)=h(nI)=h(0), h(5), h(10), h(15), h(20), h(25), h(30), h(35), h(40), h(45)p1(n)=h(1+nI)=h(1), h(6), h(11), h(16), h(21), h(26), h(31), h(36), h(41), h(46)p2(n)=h(2+nI)=h(2), h(7), h(12), h(17), h(22), h(27), h(32), h(37), h(42), h(47)p3(n)=h(3+nI)=h(3), h(8), h(13), h(18), h(23), h(28), h(33), h(38), h(43), h(48)p4(n)=h(4+nI)=h(4), h(9), h(14), h(19), h(24), h(29), h(34), h(39), h(44), h(49)第8章 多采样率数字信号处理 8.7 采样率转换器的采样率转换器的MATLAB实现实现MATLAB信号处理工具箱提供的采样率转换函数有upfirdn, interp, decimate, resample， 其功能简述如下。 Y=upfirdn(X, H, I, D)： 先对输入信号向量X进行I倍零值内插， 再用H提供的FIR数字滤波器(FIRDF)对内插结果滤波， 其中H为FIR DF的单位脉冲向量， FIR DF采用高效的多相实现结构。 最后按因子D抽取得到输出信号向量Y。 第8章 多采样率数字信号处理 Y=interp(X, I)： 采用低通滤波插值法实现对序列向量X的I倍插值， 其中的插值滤波器让原序列无失真通过， 并在X的两个相邻样值之间按照最小均方误差准则插入I1个序列值。 得到的输出信号向量Y的长度为X长度的I倍。 Y=decimate(X, D, N)： 先对序列X抗混叠滤波， 再按整数因子D对序列X抽取。 输出序列Y的长度是X长度的1/ D。 抗混叠滤波用N阶切比雪夫型低通滤波器， 阻带截止频率为s/(2D)， 如果省略N， 则默认用8阶切比雪夫型低通滤波器； Y = decimate(X, D, N, FIR)用长度为N的FIR滤波器， FIR滤波器是抽取函数decimate自动调用fir1(N, 1/D)设计的。 省略N， 则默认用30点FIR数字滤波器。 第8章 多采样率数字信号处理 Y=resample(X, I, D)： 采用多相滤波器结构实现按有理数因子I/ D的采样率转换。 如果原序列向量X的采样频率为Fx， 长度为Lx， 则序列Y的采样频率为Fy=(I/D)Fx ， 长度为(I/D)Lx（当(I/D)Lx不是整数时， Y的长度取不小于(I/D)Lx的最小整数）。 该函数具有默认的抗混叠滤波器设计功能， 按照最小均方误差准则调用函数firls设计。 Y, B =resample(X, I, D)： 返回输出信号向量Y和抗混叠滤波器的单位脉冲序列向量B。 Y=resample (X, I, D, B)： 允许用户提供抗混叠滤波器的单位脉冲响应序列向量B。 这些函数的其他调用格式请用help命令查阅。 第8章 多采样率数字信号处理 【例例】 编写程序产生长度为41的序列x(nnn)， 再调用resample函数对x(n)按因子3/8进行采样率变换， 并绘制采样率变换器的输入序列x(n)、 输出序列y(n)和采样率变换器中的FIR数字滤波器的单位脉冲响应h(n)及其频率响应特性曲线。 解解 本例题的实现程序如下： 第8章 多采样率数字信号处理 例8.7.1 实现程序ep871.m: 调用resample函数实现按因子3/8进行采样率变换n=0:40; xn=sin(0.1*pi*n)+0.5*sin(0.5*pi*n); %产生长度为41的序列向量xnyn, hn=resample(xn, 3, 8); %对xn按因子3/8进行采样率变换, yn为转换器输出序列, hn是FIRDF的单位脉冲响应%以下是绘图部分subplot(3, 2, 1); stem(hn, .); axis(0, 160, ); title(a); xlabel(i); ylabel(h(i)w=(0:1023)*2/1024; 第8章 多采样率数字信号处理 subplot(3, 2, 2); plot(w, 20*log10(abs(fft(hn, 1024); axis(0, 1/4, 80, 20); grid ontitle(b); xlabel(omega/pi); ylabel(20lg(|Hg(omega)|)subplot(3, 1, 2); stem(n, xn, .); title(c); xlabel(n); ylabel(x(n)ny=0:length(yn)1; subplot(3, 1, 3); stem(ny, yn, .); title(d); xlabel(m); ylabel(y(m)ZK) 第8章 多采样率数字信号处理 运行程序得到采样率转换器的输入信号x(n)、 输出信号y(m)、 FIR数字滤波器的单位脉冲响应h(n)及其频率响应的波形分别如图8.7.1（c）、 （d）、 （a）和（b）所示。 由图（d）可以看出， resample函数默认设计的抗混叠滤波器的阻带截止频率为/8， 阻带最小衰减大于40 dB, 满足理论要求。 x(n)的长度为40， y(m)的长度为403/8=15， h(n)的长度为161。 由图8.7.1可见， 低通滤波器滤除输入信号中较高的频率成分0.5sin(0.5n)， 采样率变换器输出单频正弦波采样， 采样频率为Fy=3Fx/8， 所以输出序列第8章 多采样率数字信号处理 图8.7.1 程序ep871.m运行结果第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. 假设信号x(n)及其频谱X(ej)如题2图所示。 按因子D=2直接对x(n)抽取， 得到信号y(m)=x(2m)。 画出y(m)的频谱函数曲线， 说明抽取过程中是否丢失了信息。 第8章 多采样率数字信号处理 题2图第8章 多采样率数字信号处理 3. 按整数因子D=4抽取器原理方框图如题3图（a）所示。 其中， Fx=1 kHz, Fy=250 Hz, 输入序列x(n)的频谱如题3图（b）所示。 请画出题3图（a）中理想低通滤波器 hD(n)的频率响应特性曲线和序列v(n)、 y(m)的频谱特性曲线。 4. 按整数因子I内插器原理方框图如题4图所示。 图中， Fx=200 Hz, Fy=1 kHz， 输入序列x(n)的频谱如题3图（b）所示。 确定内插因子I， 并画出题4图中理想低通滤波器hI(n)的频率响应特性曲线和序列v(m)、 y(m)频谱特性曲线。 第8章 多采样率数字信号处理 题3图第8章 多采样率数字信号处理 题4图第8章 多采样率数字信号处理 5*. 设计一个抽取器， 要求抽取因子D=5。 用remez函数设计抗混叠FIR滤波器， 图示滤波器的单位脉冲响应和损耗函数。 要求通带最大衰减为0.1 dB， 阻带最小衰减为30 dB， 过渡带宽度为。 画出实现抽取器的多相结构， 并求出多相实现时各子滤波器的单位脉冲响应。 6*. 设计一个内插器， 要求内插因子I=2。 用remez函数设计镜像FIR滤波器， 图示滤波器的单位脉冲响应和损耗函数。 要求通带最大衰减为0.1 dB， 阻带最小衰减为30 dB， 过渡带宽度为。 画出实现内插器的多相结构， 并求出多相实现时各子滤波器的单位脉冲响应。 第8章 多采样率数字信号处理 7*. 设计一个按因子2/5降低采样率的采样率转换器， 画出系统原理方框图。 要求其中的FIR低通滤波器通带最大衰减为1 dB， 阻带最小衰减为30 dB。 设计FIR低通滤波器的单位脉冲响应， 并画出一种高效实现结构。8*. 假设信号x(n)是以奈奎斯特采样频率对模拟信号 xa(t)的采样序列， 采样频率Fx=10 kHz。 现在为了减少数据量， 只保留0f3 kHz的低频信息， 希望尽可能降低采样频率， 请设计采样率转换器。 要求经过采样率转换器后， 在频带0f2.8 kHz中频谱失真不大于1 dB， 频谱混叠不超过1。 第8章 多采样率数字信号处理 （1） 确定满足要求的最低采样频率Fy和相应的采样率转换因子。 （2） 画出采样率转换器原理方框图。 （3） 确定采样率转换器中FIR低通滤波器的技术指标， 用等波纹最佳逼近法设计FIR低通滤波器， 画出滤波器的单位脉冲响应及其损耗函数曲线， 并标出指标参数（通带截止频率、 阻带截止频率、 通带最大衰减和阻带最小衰减）。 （4） 求出多相实现结构中子滤波器的单位脉冲响应， 并列表显示或打印。