1. 多速率信号处理动手仿真前先重新认识采样率转换先说结论多速率信号处理不是把采样率改一改这么简单它牵扯到频谱结构的变化、滤波器的边界设计以及实时系统里那点算力到底够不够花。很多初学者拿到滤波器设计与仿真_13.多速率信号处理这个题目第一反应是去翻MATLAB的resample函数跑通一条曲线就觉得完事了——这恰恰是最大的坑。多速率信号处理解决的核心问题只有一个在数字域内改变信号的采样率并且让改变之后的信号在频谱意义上“干净、无损”。它广泛出现在软件无线电、音频重采样、声呐处理、图像缩放和通信基带处理里。举个例子你手机里把48kHz采样率的音频转成44.1kHz输出给声卡这就是多速率处理你在FPGA里把ADC采回来的100MHz数据降成25MHz再进基带算法这也是多速率处理。前者叫分数倍速率转换后者叫整数倍抽取名字不同底层逻辑一致。这篇内容适合三类人正在做信号处理课程设计的在校生刚接触软件无线电或音频DSP的工程师以及所有想彻底搞懂“为什么降采样之前必须低通滤波”这个问题的人。我默认你具备最基础的信号与系统知识比如傅里叶变换、抽样定理和FIR滤波器的基本概念。如果你连这些也还没完全消化我会在讲到频谱变化时尽量用直观的方式带一遍。我先用一句话把整个章节的思路串起来多速率信号处理 抽取/插值 配套滤波器 高效实现结构。抽插决定“变多少”滤波器决定“质量好不好”多相结构决定“能不能跑起来”。这三件事分开理解都不难难的是把它们组合到同一个工程场景里去。下面我把每一步拆开讲重点说仿真参数怎么定、滤波器怎么选、实现结构怎么搭以及我在实际调试中踩过的那些坑。2. 核心概念与频谱视角抽取、插值究竟动了什么2.1 抽取降采样的本质是“频谱复制”不是“压缩”假设原始信号采样率是 fs抽取因子为 M也就是每隔 M-1 个点丢掉 M-1 个保留 1 个输出采样率变成 fs/M。这个操作在时域上很好理解但频域上发生的事情很容易被忽略。在数字域信号的频谱是周期的周期为 fs。抽取之后的新采样率是 fs/M新的周期也是 fs/M。如果原信号在 (fs/(2M), fs/2) 区间内存在非零频谱分量这些分量在抽取后就会被“折叠”回 (0, fs/(2M)) 区间内和原本低频段的信号混在一起产生无法去除的混叠。我在给学生演示这个现象时最喜欢用两个正弦叠加的例子一个 100Hz一个 3200Hz采样率 8000Hz。直接做 4 倍抽取后3200Hz 分量会混叠到 800Hz 附近——因为 3200 对 2000即 8000/4求模后等于 200再加上原始折叠关系实际出现在 800Hz。这根本不是你想要的结果。所以在抽取操作之前必须用一个低通滤波器把截止频率限制在 fs/(2M) 以内。这个滤波器通常叫抗混叠滤波器。仿真时很多人直接调用downsample或resample内部会自动加滤波器但你自己设计系统时滤波器的阶数、通带纹波、阻带衰减都要明确。2.2 插值升采样的“镜像”问题插值是在每两个原始采样点之间插入 L-1 个零值采样率变为 L×fs。插零之后频谱中除了原始信号频带0, fs/2还会在 (L×fs/2 附近) 出现原始频谱的周期复制镜像。这些镜像必须用低通滤波器滤除否则输出信号听起来会有“金属声”信号分析时会看到一堆不该存在的高频分量。这个滤波器叫抗镜像滤波器。关键点在于抽取前要限带插值后要滤镜像两个滤波器看似方向相反本质上都是低通滤波截止频率取两者中较小的那个。如果是分数倍变换 L/M滤波器的归一化截止频率应该取 min(pi/L, pi/M)。2.3 为什么不能直接改采样率当成重采样这里要特别强调一个致命误区很多学习者在搞不太清楚的情况下直接用线性插值或者interp1去“重采样”信号。这在时域上能出来一条曲线但频谱上完全不可控。线性插值本质上是一个很差的低通滤波器它的频响在接近奈奎斯特频率时并不理想会带来衰减和失真。对于音频听感测试可能勉强能忍但做通信信号分析、雷达信号处理这种方案直接废掉。我可以给一个直观对比表方案可处理的场景风险复杂度直接抽取/插值仅理论演示不做真实信号处理严重混叠/镜像低线性插值重采样简单波形显示高频严重失真低FIR抽取/插值绝大多数高质量处理场景设计好滤波器则无中多相结构实时系统、资源受限场景设计稍复杂中高CIC补偿滤波极高倍率变化的前后端需补偿通带中表格最后一行的CIC滤波器是专门应对极大倍数抽取/插值场景的后面我会重点展开。现在先记住一个原则凡是说“重采样”的地方心里默认它包含了滤波器设计而不是简单的插值/抽点。3. 滤波器设计与参数计算绕不开的抗混叠/抗镜像低通3.1 设计前的参数映射在设计滤波器之前先把几个关键参数定义清楚。我以工程里最常见的场景来举例系统输入采样率 fs_in 48kHz需要转换成 fs_out 16kHz也就是 3 倍抽取。需要设计一个低通滤波器在抽取前使用。三个设计参数最核心通带截止频率 fp决定了保留信号的有效带宽。假设你只关心音频里 5kHz 以下的部分那 fp 可以设为 5kHz通带纹波对应人耳可感知的幅度波动一般取 0.1dB 以内。阻带起始频率 fst决定了哪里开始抑制。为了不混叠必须保证在折叠频率 fs_out/2 8kHz 以上完全抑制。也就是说 fst 至少要取 8kHz且 8kHz 以上的衰减要足够大。阻带衰减 As工程上通常取 60dB这个值是绝大多数信号处理系统的底线。低于 40dB 会在抽取后留下可见的混叠噪声。根据这三个值用MATLAB的firpm或fdesign.lowpass来设计。一个比较省钱的思路是直接用fdesign对象自动算出所需的最小阶数。我在实际项目里从来不会手动去选filterOrder——让设计工具自动算阶数比人眼估靠谱得多。3.2 计算过程实例48kHz转16kHz的抗混叠滤波我带着大家一起算一遍。已知fs_in48000M3信号有效带宽 0~5000Hz阻带从 8000Hz 开始阻带衰减要求 80dB给自己留余量。归一化频率计算通带边界fp / (fs_in/2) 5000 / 24000 ≈ 0.2083归一化到奈奎斯特频率阻带起始fst / (fs_in/2) 8000 / 24000 ≈ 0.3333设计代码fsIn 48000; M 3; fp 5000; fst 8000; As 80; d fdesign.lowpass(Fp,Fst,Ap,As, fp, fst, 0.1, As, fsIn); h design(d, equiripple); disp(length(h.Numerator)); % 查看滤波器阶数实际设计出来可能是一个两百多阶的FIR。两百阶对于48kHz采样率、3倍抽取来说完全够用。为什么FIR在这里是首选而不是IIR因为多速率处理里滤波器必须保持严格线性相位IIR虽然阶数低但相位是非线性的在通信和音频场景里会引发群延迟失真。这一点几乎可以视为铁律。3.3 滤波器设计质量验证仿真前必须看的三张图设计完滤波器先别急着接进系统里。用freqz看幅频响应同时做三件事检查通带纹波是否真的在 0.1dB 范围内——买要紧的窄带信号纹波太大会造成幅度失真。检查阻带衰减是否达标——用yline(-80)画一条参考线阻带响应整体压在这条线以下。检查过渡带宽度是否满足 5000Hz 到 8000Hz 的要求——如果过渡带太宽阻带起始频率被迫提前就会牺牲有效带宽。我见过太多人设计完滤波器直接fvtool(h)瞄一眼就跳到下一步最后发现输出信号里带着残余混叠噪声却找不到原因。把这三项验证固化到自己的调试流程里能省掉后面一大堆排查时间。3.4 级联抽取与多级实现的算力对比当抽取倍数 M 很大比如从 96kHz 直接降到 8kHz也就是 12 倍抽取如果一步到位设计滤波器阶数会非常高实时系统里乘法器资源直接爆炸。解决方案是多级抽取先降 2 倍滤波再降 2 倍滤波再降 3 倍滤波。每一步的滤波器阶数都显著低于单级设计整体乘法计算量大幅下降。我拿一个实际例子对比单级 48 倍抽取需要大约 1000 阶 FIR拆成 2、3、4、2 四级抽取后每一级滤波器大约 50~200 阶总计算量减少了约 70%。这就是“多速率”中“多”的另一层含义不只是一次变速率而是把速率变换拆成多级用算力换性能、用结构换资源。4. 仿真与实现从原型滤波器到多相结构再到CIC4.1 第一步先搭一个不含滤波器的“裸抽取”仿真找不到好例子去理解过滤跟不过滤的区别最好的办法就是亲自跑一把。生成一个含 1kHz 和 7kHz 两个正弦分量的信号采样率 16kHz做 2 倍抽取。fs 16000; t (0:999)/fs; x sin(2*pi*1000*t) 0.5*sin(2*pi*7000*t); M 2; ydown downsample(x, M); f_ydown linspace(0, fs/2, length(ydown)); Y fft(ydown); plot(f_ydown, 20*log10(abs(Y(1:length(ydown)/2))));未经过滤波器的话7kHz 分量会混叠到 1kHz 附近从频谱上会看到一个“不明真相的突起”。这就是混叠噪声。做仿真的价值就在于此你可以“眼见为实”地看到理论公式在说什么。4.2 第二步加入滤波器跑通完整抽取链路用上一章设计的滤波器 h执行完整流程xf filter(h, x); % 抗混叠滤波 y downsample(xf, M); % 抽取这里有个细节filter和conv的结果在稳态后相同但起始段会因为滤波器延迟而产生差别。实际工程中滤波器群延迟会导致输出信号相对原始信号有固定的时间偏移这个偏移等于滤波器阶数的一半。如果你在意时序对齐需要做补偿。4.3 第三步插值链路与插值后滤波插值侧和抽取侧对称但先后顺序不同先插零再滤波。L 2; xu upsample(x, L); % 插零 y filter(h, xu); % 抗镜像滤波多相结构在实现上的意义就在于它把插零后那大量无意义的零乘运算直接跳过了只在有效位置做乘法。这一步对实时系统很关键。4.4 多相分解的 MATLAB 实现多相分解是高效实现的核心。以 3 倍抽取为例把原型滤波器 h 按相位拆成 3 个子滤波器。在 MATLAB 里可以用reshape实现M 3; hphases reshape(h, M, []).; % 每列是一个子滤波器可能长度不等实际使用时如果滤波器长度不是 M 的整数倍需要补零到整数倍否则 reshape 会报错。这是我第一次写多相滤波时踩过的坑印象极深。补零操作会产生相位偏移所以补零量要计入群延迟补偿。对多速率系统做MATLAB仿真时可以用 dsp 工具箱里的dsp.FIRInterpolator和dsp.FIRDecimator它们内部已经用多相结构实现你不必手写。但如果要移植到FPGA或者嵌入式平台手动实现多相结构就是必修课。4.5 CIC滤波器极大量率变化时的“暴力”方案CIC级联积分梳状滤波器是另一种思路不含乘法器只做加减法在FPGA里极其节省资源。它特别适合做高倍率抽取的前端比如先把 100MHz 降到 25MHz再做精细滤波把 25MHz 降到 100kHz。CIC的传递函数等效于一系列积分器和梳状器的级联。它的频率响应在通带内有明显下垂所以通常需要再接一个FIR补偿滤波器把通带拉平。MATLAB里用fdesign.cic直接设计有补偿的CICcicobj fdesign.cic(3, 16); % 3级CIC抽取因子16 f design(cicobj);CIC 设计有一个关键选择级联级数 N。N 越大阻带衰减越好但通带下垂越严重。工程上 N3~5 是常见取值。我在一个声呐项目里用过 4 级 CIC配合一个 128 阶补偿 FIR整体效果比纯 FIR 方案节省了接近 90% 的乘法器资源。5. 仿真常见问题与排查实录5.1 问题一抽取后频谱里出现了“不该有”的低频分量原因通常是抗混叠滤波器阻带衰减不足或者截止频率设置过高。排查方法计算混叠频率位置看该处是否确实存在峰。具体步骤是先在滤波器设计阶段确保阻带衰减在目标频段达到 60dB 以上。如果滤波器阶数受限考虑增大过渡带宽度。5.2 问题二插值后信号高频端有明显衰减这通常是抗镜像滤波器的过渡带太宽导致镜像频率未被完全抑制反而在目标频段留下了衰减痕迹。我遇到过一次用了一个通带截至到 0.45×Nyquist 的滤波器做过 L2 插值结果输出 15kHz 附近衰减了 2dB。检查后发现滤波器阻带起始频率设太高了镜像频率落在阻带边缘。后来把阻带起始降到 0.47×Nyquist衰减恢复了。5.3 问题三实时实现时乘法器资源爆表第一个思路是换用多相结构只在高有效相位上做运算。第二个思路是级联抽取把一次大变化拆成多次小变化。第三个思路是评估 CIC。在FPGA上我一般是 CIC 做粗降FIR 做精补这样资源最省。用 MATLAB 的filtersupport类可以估算滤波器实现成本不过我在实际项目里更多的是直接看综合报告里的 DSP48 使用数量。5.4 问题四分数倍变换L/M的序列顺序搞反正确顺序是先插 L 倍再滤波再抽取 M 倍。如果先抽取后插值会丢失信息后续插值只能“脑补”缺失的高频分量。这一点是顺序错误的高频雷区务必记住。5.5 问题五群延迟导致输出信号时域错位任何 FIR 滤波器都有固定群延迟等于 (N-1)/2 个采样点。在多速率系统中抽取或插值会改变这个延迟对应的输出采样位置。如果你做的是需要严格时基对齐的测量类应用比如测距雷达或声学回声消除必须把群延迟补偿纳入系统总延时预算。我用过一个简单方案把滤波器用grpdelay函数计算整数群延迟然后在输出端裁剪固定的采样点。这个方法在整数倍速率转换时非常稳定。5.6 问题六定点仿真浮点仿真结果不一样这是往硬件移植时逃不掉的问题。浮点仿真能跑通不代表定点也能过尤其是滤波器系数量化后阻带衰减可能下降 10~20dB。正确的做法是先用fipref和fi对象做定点建模观察量化前后的频谱差异。另外一个呼声很高的问题是CIC 溢出了怎么办CIC 的积分器是无条件增长的需要根据输入位宽和级联级数扩展寄存器位宽这也是定点设计中必须算清楚的一件事。6. 心得与补充把多速率仿真变成工程直觉多速率信号处理这个章节在学习路径上的地位有点像“分水岭”。在这之前你处理的都是单个采样率下的线性系统在这之后你开始在多个速率域之间穿梭采样率不再是固定的、理所当然的前提条件。这一步跨过去后续理解软件无线电架构、数字下变频、欠采样、滤波器组时都会顺很多。我个人在带项目时发现能把多速率仿真真正做明白的同学通常具备两个共同点一是愿意从频谱图上去解读每一步操作的后果不把resample当黑盒二是愿意动手写多相结构的代码哪怕只是仿真层面也要把“零乘被跳掉”这个行为亲手验证一遍。一些实操建议按性价比排序建议一在仿真阶段始终保留滤波前后的频谱对比图。不要只画最终输出。多速率系统里错误发生在哪一级往往只有频谱能告诉你。建议二滤波器阶数自动设计出来很高别急着砍先用多级方案替代。在很多情况下多级滤波器总阶数比单级少一个数量级以上。建议三做分数倍变换时先用最接近的整数倍去估算性能再调整到精确比值。这样调试起来心智负担小很多。建议四仿真环境里把采样率和滤波器的归一化频率对应关系写清楚我习惯在代码头部输出一行校验信息确保设计参数没有搞错量级。另外一个经常被忽略的细节是多速率处理后输出信号的功率会发生变化。抽取会降低噪声功率但也会改变信号功率水平插值同样如此。所以做信噪比分析时必须按采样率变化修正功率比例否则算出来的 SNR 有系统性偏差。如果你做的是音频重采样还要注意听感测试比频谱测试更重要。频谱干净不代表听感好尤其在处理瞬态信号时滤波器群延迟的色散特性会造成类似“回声”的听感。此时可以尝试线性相位 FIR 和最小相位滤波器进行 AB 对比。再补充一个我在项目里常用的工具链MATLAB 负责算法验证和滤波器设计Python 的scipy.signal.resample_poly负责快速原型对比。两边结果一致后我再把滤波器系数导出成查找表进入 FPGA 定点验证阶段。这个流程我已经稳定用了好几年换个角度说它也保障了移植阶段最大程度的确定性。多速率仿真并不是一个孤立的技术点它是整个数字信号处理系统中连接“前端采样”和“后端算法”的桥梁。你把这个桥梁设计得越扎实后续的信号质量就越好。这一章虽然排在第十三节但它背后的思想会贯穿到后面所有章节里。希望这篇拆解能帮你在做仿真时少走几条弯路。如果哪一步的计算或设计还有疑问最好的办法是打开 MATLAB把代码跑一遍亲眼看看频谱上发生的变化。多速率处理这个主题靠读是读不会的亲手改一个参数、观察一次频谱跳动比背十遍公式都管用。
