数据压缩课最麻烦的作业往往不是写熵编码或者量化器而是让你把一整套滤波器组从原理到仿真完整实现出来。完全重建QAM滤波器组这个题目表面上是调库调用实际上牵扯到多速率信号处理、原型滤波器设计、调制映射、延迟对齐和混叠分析每个环节都能埋雷。这篇文章把我做这个作业时踩过的坑、总结过的推导、以及可以直接跑通的实现步骤都整理出来给正在被这门作业折磨的人一个参考。1. 作业到底在做什么QAM滤波器组与子带编码1.1 数据压缩作业里的“滤波器组”角色数据压缩的核心思路无非是去掉冗余。时域预测编码去的是时间冗余变换编码去的是空间或统计冗余而滤波器组做的事情是帮你在频域上把信号切成若干个子带然后根据每个子带的能量和感知重要程度分配合适的量化位数。这个过程在语音编码、音频编码和图像压缩里都有应用经典如子带图像编码、MDCT音频编码背后都是滤波器组的思路。你可能会问既然有FFT这种“一刀切”的变换为什么还要用滤波器组区别在于FFT是对整段信号做固定分块块与块之间没有任何平滑过渡会导致很明显的边界效应。而滤波器组的设计天然带有重叠和过渡分析端把信号划分到子带综合端又能把子带重新拼回完整信号中间的切换是连续的重构质量就靠“完全重建”这个性质来保证。作业里提到的QAM滤波器组就是其中一种实现方式。它名字里带着“QAM”容易让人联想到通信里的正交幅度调制实际上它跟通信里的星座图没有直接关系。这里的QAM是指用复指数调制正交调制把原型低通滤波器搬移到不同频带从而形成一组相互覆盖整个频带的子带滤波器。这种设计的好处是结构统一分析端和综合端的滤波器关系明确配合上过采样处理可以做到工程意义上的完全重建。1.2 QAM滤波器组的名词拆解与整体结构先把滤波器的组立结构讲清楚。一个标准的M通道滤波器组包含分析滤波器和综合滤波器两部分。分析端有M个滤波器编号从0到M-1。每个分析滤波器 $H_k(z)$ 的作用是从原始信号中提取第k个子带的成分。原始信号经过这M个滤波器后得到M路信号每一路都只保留了一个频带内的内容然后经过抽取降低采样率。综合端也有M个滤波器编号对应。M路子带信号先经过插值恢复到原始采样率再通过综合滤波器 $G_k(z)$ 滤掉插值产生的镜像谱最后把所有子带输出加起来得到重建信号。QAM滤波器组的结构特点是所有分析滤波器都由同一个低通原型滤波器 $h(n)$ 调制而来。具体做法是$$h_k(n) h(n) \cdot e^{j \omega_k (n - \tfrac{N-1}{2})}$$其中 $\omega_k$ 是第k个子带的中心频率$N$ 是原型滤波器长度。这个“调制”操作本质就是把低通滤波器的频谱搬到 $\omega_k$ 附近。$\omega_k$ 的取法有讲究好的取值能让子带之间刚好无缝衔接也方便后续综合端重建时消除混叠。用这种“一个原型调制”的方式你只需要设计好一个低通滤波器剩下的M个带通滤波器全部由公式生成调试起来很舒服。后续如果要改通道数也不用推翻重来改个M重新生成一遍调制系数就行。这个特性写作业和做工程都很省心。1.3 为什么“完全重建”是硬指标滤波做得好不好最终要落到一个验收标准上原始信号 $x(n)$ 经过分析-综合系统后输出 $\hat{x}(n)$ 与 $x(n)$ 之间差多少。如果系统是理想的那么除了固定延迟和幅度缩放输出应该和输入一模一样。这种“输入输出只差一个延迟和常数增益”的性质就叫完全重建。对数据压缩而言完全重建意味着只要量化环节不做任何有损操作或者量化误差为零整个编解码链路就不会引入额外失真。在这个前提下压缩的误差完全由量化步长决定系统设计师可以放心地调整压缩参数而不需要担心滤波器组本身在偷偷制造噪声。但现实是滤波器组涉及抽取和插值必然会产生混叠。如果原型滤波器设计得不好或者调制参数不对混叠就会泄漏到重构信号中表现出来就是高频噪声或者周期性的波纹。这些问题你在时域里看波形可能不明显但在频域里一看就露馅了。所以判断作业做得对不对先检查重构误差再检查混叠泄漏。这两个指标合格才算真正把完全重建做出来了。2. 完全重建背后的三个失真项2.1 分析-综合系统的输入输出关系要理解完全重建先要把整个链路的数学表达式写出来。这里我们假设抽取因子为 $D$插值因子也是 $D$。输入端信号 $X(z)$ 经过第k个分析滤波器 $H_k(z)$ 后被 $D$ 倍抽取再被 $D$ 倍插值最后经过综合滤波器 $G_k(z)$ 后求和。这个过程中抽取会产生频谱折叠。只要分析滤波器的阻带不是零点抽取后的子带信号中就会混入相邻频带的成分。插值会产生镜像谱需要在综合端通过滤波器抑制。把这些效应都考虑进去后输出可以表示为$$Y(z) \frac{1}{D}\sum_{l0}^{D-1} X(zW_D^l) \sum_{k0}^{M-1} G_k(z) H_k(zW_D^l)$$其中 $W_D e^{-j2\pi/D}$。这个式子看着吓人其实拆开看就两件事求和符号里 $l0$ 的那一项对应的是原始信号的正常传输路径而 $l \neq 0$ 的各项对应的是因为抽取折叠而产生的混叠成分。每一项混叠的强度取决于分析滤波器在折叠频点处的响应和综合滤波器在该处的抑制能力。为了简化表达定义传输函数$$T_l(z) \frac{1}{D}\sum_{k0}^{M-1} G_k(z) H_k(zW_D^l)$$那么输出就是$$Y(z) T_0(z) X(z) \sum_{l1}^{D-1} T_l(z) X(zW_D^l)$$到这里完全重建的三个条件就浮出水面了。2.2 混叠、幅度失真、相位失真分别如何产生混叠失真对应的是 $T_l(z) 0$对所有 $l \neq 0$ 成立。也就是说那些因为抽取而折叠进来的频谱成分必须在综合端被彻底抵消。理想情况下分析滤波器的阻带足够干净抽取后折叠进来的成分本来就很小但更重要的一点是综合滤波器的频响覆盖范围要能补回折叠丢失的部分。在QAM滤波器组中由于各子带滤波器是同一个原型调制得到的相邻子带的尾部响应可以相互抵消这正是调制结构比任意滤波器组的优势所在。幅度失真对应的是 $T_0(z)$ 的幅频响应必须是常数或者近似常数。如果 $T_0(z)$ 在某些频率处有明显的谷或峰那么原信号经过系统后某些频段会被削弱或放大听起来就是“频谱不平坦”。对于完全重建系统要求 $|T_0(e^{j\omega})|$ 在所有频率上保持一致通常归一化为1。相位失真对应的是 $T_0(z)$ 的相位必须与频率成线性关系也就是整体只引入一个固定延迟 $\tau$。如果相位呈现非线性信号的不同频率分量会被错开波形就会产生畸变。线性相位原型滤波器在这里帮了大忙因为只要原型 $h(n)$ 是线性相位的调制后的分析滤波器也保持线性相位特性综合端再取合适配对整体链路的相位就可控。这三个条件都满足时$Y(z) c z^{-\tau} X(z)$也就是完全重建。实际操作中工程意义上的完全重建允许极小的偏差比如重构误差低于-60dB就可以当作“重建成功”来对待。2.3 过采样QAM为何能实现完全重建你可能在教科书里看到过临界采样的M通道滤波器组要做完全重建原型滤波器的设计条件非常苛刻尤其是M比较大的时候几乎很难同时压住混叠和幅度失真。那为什么QAM滤波器组还能做关键在于它可以配合过采样使用。过采样的含义是分析滤波器个数是M但抽取因子D小于M。最常见的配置是 $D M/2$也就是2倍过采样。原本信号被切成M个子带每个子带只保留 $1/M$ 的带宽信息按理说采样率可以降为原来的 $1/M$但过采样只降到 $1/D$等于每个子带信号里留了余量频谱之间留有保护带抽取折叠时混叠就小得多。打个比方临界采样就像把一整条面包切成M片每片都紧挨着切歪一点就露馅。过采样则是让每片之间留一点空隙切成M片但每片都薄一点拼接的时候自然不容易出错。对于完全重建来说保护带就是容错空间。QAM调制再加上过采样综合端重建时那些 $l \neq 0$ 的混叠项会被原型滤波器的阻带衰减和过渡带设计联合压制重构质量可以做得非常高。这也是为什么作业里推荐使用过采样QAM而不是临界采样QMF前者在工程实现上简单得多对原型滤波器的要求也更宽松。理解了这点后面设计参数时你就知道哪些地方可以放松哪些地方必须死磕。3. 原型滤波器设计与调制构造步骤3.1 设计M通道系统的关键参数做QAM滤波器组第一步是定参数。参数之间互相制约定错了后面返工成本很高。我整理了一张参数表写作业和做毕设都可以直接参考。参数符号取值依据我的推荐通道数M子带数量取决于压缩需求M8或M16量级适中抽取因子D过采样率决定 DMDM/2典型2倍过采样原型滤波器长度N越长过渡带越窄但延迟越大N4KM-1K取2~4原型截止频率fc要覆盖子带占用带宽fc1/(2*D)归一化频率阻带衰减AS决定混叠泄漏水平60dB以上最好80dB过渡带宽Δf与N直接相关由AS和N反推这里尤其要注意原型滤波器长度和过渡带宽的关系。窗函数法设计滤波器时过渡带宽与滤波器长度近似成反比。N太小过渡带宽就宽相邻子带会互相污染N太大整个系统的延迟也大而且接近长滤波器的数值精度问题会浮现出来。我的经验是先定M再定D然后用 $N 4KM - 1$ 作为初始值如果重构误差不够再把K往上调。截止频率的取值也别拍脑袋。分析滤波器的通带要覆盖 $\pi/M$ 范围内的信息而抽取因子D决定了子带信号保留的带宽。当 $D M/2$ 时每个子带的保留带宽约为 $1/(2D) 1/M$所以原型滤波器的通带截止频率取 $\pi/M$ 附近比较合理。你可以把原型滤波器的6dB衰减点设置在这个位置保证相邻子带的交叠处能量互补。3.2 用窗函数法设计原型低通滤波器确定好参数后设计原型滤波器最省事的方法就是窗函数法尤其是Kaiser窗因为它允许你独立控制阻带衰减和过渡带宽不像矩形窗那样参数是死的。Kaiser窗设计低通滤波器的步骤如下确定阻带衰减 $A_s$单位dB比如取80dB。根据衰减值计算Kaiser窗的 $\beta$ 参数$$\beta \begin{cases} 0.1102(A_s - 8.7), A_s 50 \ 0.5842(A_s - 21)^{0.4} 0.07886(A_s - 21), 21 \le A_s \le 50 \end{cases}$$确定过渡带宽 $\Delta\omega$然后反推滤波器长度$$N \frac{A_s - 7.95}{2.285 \cdot \Delta\omega} 1$$过渡带宽要根据子带间距来定。子带之间能容忍的过渡带范围大约是 $0.1\pi/M$ 到 $0.4\pi/M$再宽就会侵占相邻子带的频带。用MATLAB做这一步非常直接。下面是我作业里实际用过的设计代码function h design_prototype(M, D, K, As) % 设计Kaiser窗原型低通滤波器 % M: 通道数 % D: 抽取因子 % K: 长度系数滤波器长度为 4*K*M - 1 % As: 阻带衰减 (dB) N 4 * K * M - 1; fc 1 / (2 * D); % 由阻带衰减计算Kaiser窗beta参数 if As 50 beta 0.1102 * (As - 8.7); else beta 0.5842 * (As - 21)^0.4 0.07886 * (As - 21); end % 理想低通脉冲响应线性相位 n 0:N-1; h_ideal 2 * fc * sinc(2 * fc * (n - (N-1)/2)); % 加Kaiser窗 w kaiser(N, beta); h h_ideal .* w; % 归一化使各子带合成增益为常数 h h / sum(h); end这里有一个细节很多教程不提为什么长度取 $4KM-1$ 而不是随便一个偶数因为QAM滤波器组在做调制时需要原型滤波器关于中心点对称即满足 $h(n) h(N-1-n)$这要求N是奇数。同时长度要是M的倍数关系这样调制后的相位关系才规整综合端重建时各子带之间的相位对齐才有保证。$4KM-1$ 这个形式保证了N既是奇数又与M有明确的模关系后续分析延迟和做多相分解都会省事很多。3.3 从原型到QAM滤波器组调制这一步怎么算原型滤波器设计好后剩下的就是调制。QAM滤波器组使用复指数调制第k个子带的中心频率为$$\omega_k \frac{2\pi(k 0.5)}{M}, \quad k 0, 1, \dots, M-1$$这里加了一个 $0.5$ 的偏移不是从0开始而是从半整数倍开始。这样做的目的是让第一个子带的中心频率离开0频避免DC分量对第一个子带产生干扰同时也让原型滤波器在0频附近的响应不必完美因为第一个子带不再覆盖直流点。分析滤波器的生成公式$$h_k(n) h(n) \cdot e^{j\omega_k (n - \tfrac{N-1}{2})}$$注意调制指数里的 $n - (N-1)/2$ 这个平移项。加上它是为了保证 $h_k(n)$ 仍然保持关于中心点的共轭对称性也就是$$h_k(N-1-n) h_k^*(n)$$这个性质在后面做延迟补偿时很有用。写代码的时候平移项千万别漏否则滤波器的相位会变得很乱重构出来的信号波形会歪得离谱。生成所有分析滤波器的MATLAB代码function [H, G] design_qam_filterbank(h, M) % 基于原型滤波器h生成M通道QAM分析滤波器组H和综合滤波器组G N length(h); n (0:N-1) - (N-1)/2; % 中心对齐的时间索引 H zeros(M, N); G zeros(M, N); for k 0:M-1 w_k 2 * pi * (k 0.5) / M; mod exp(1j * w_k * n); H(k1, :) h .* mod; % 综合滤波器取共轭对称形式保证相位对齐 G(k1, :) conj(fliplr(H(k1, :))); end end综合滤波器取分析滤波器的共轭时间反转这个动作后面会解释。这里先记住分析端的滤波器负责把信号拆开综合端的滤波器负责把拆开的信号还原两者是配对关系不是随意组合。3.4 分析/综合滤波器配对与延迟补偿综合滤波器取共轭时间反转是有讲究的。因为原型滤波器是线性相位的调制后的分析滤波器具备共轭对称性。综合滤波器取 $G_k(z) z^{-(N-1)} H_k^(1/z^)$在实数系下等价于对时间索引反转并取共轭。这样做可以保证整个分析-综合链路的传递函数中非延迟项相互抵消只留下一个整体延迟。如果综合滤波器随便取比如直接用分析滤波器本身而不做任何反转那么链路中每个子带的相位变化会叠加输出信号会发生频率相关的相位扭曲。在作业的验证阶段你可能发现重构波形幅度差不多但过零点位置全乱了那就是配对方式出了问题。延迟补偿也是关键。滤波器组整体链路是有通带延迟的对于长度N的原型滤波器分析滤波过程引入的延迟大约是 $(N-1)/2$综合过程同样引入 $(N-1)/2$再加上抽取插值处理带来的处理延迟 $D-1$ 等总的输出延迟大概为$$\tau \frac{N-1}{2} \frac{N-1}{2} (D-1) N - 1 D - 1$$这个公式在不同实现里会有细微差别和你滤波器的边界处理方式有关。最稳妥的办法不是手工推算而是用单位冲激做一次“冲激响应测试”直接测量输出峰值位置得到准确延迟。具体做法在第4章讲。4. 用仿真验证完全重建的完整流程4.1 搭建分析滤波器组与综合滤波器组参数设计完了接下来就是把整个系统跑起来。我建议用一个测试脚本串联整个仿真生成输入信号、通过分析端、过采样处理、通过综合端、对比重构误差。这样改任何参数跑一遍就能看到结果。分析端的处理流程是输入信号 $x(n)$ 分别与M个分析滤波器做卷积得到M路滤波结果然后每路按抽取因子D进行抽取。注意这里抽取要从正确的样本位置开始。MATLAB的downsample函数默认从第1个样本开始但滤波有延迟如果你从第1个样本就抽取等于把最前面的暂态也算进去了重构会不准。正确做法是先滤波再把前 $(N-1)/2$ 个样本丢弃或做补偿然后再抽取。综合端的流程是每路子带信号先按D倍插值补零然后与对应的综合滤波器做卷积最后把M路结果相加。function y analysis_synthesis(x, H, G, D) % 完整的分析-综合链路 % x: 输入信号行向量 % H: 分析滤波器组M x N % G: 综合滤波器组M x N % D: 抽取/插值因子 [M, N] size(H); % 计算总延迟后面用于对齐 tau_total N - 1 D - 1; % 分析端 subbands zeros(M, length(x)); for k 1:M tmp conv(x, H(k, :)); tmp tmp((N-1)/2 1 : end - (N-1)/2); % 去掉滤波暂态 subbands(k, :) tmp; end % 抽取 decimated subbands(:, 1:D:end); % 插值 upsampled zeros(M, size(decimated, 2) * D); for k 1:M upsampled(k, 1:D:end) decimated(k, :); end % 综合端 y zeros(1, size(upsampled, 2) N - 1); for k 1:M tmp conv(upsampled(k, :), G(k, :)); y y tmp; end % 截取有效部分并与输入对齐 y y(tau_total 1 : tau_total length(x)); end这个函数把抽取和插值都显式写出来了方便你观察每一步的形状。实际跑起来如果信号太长可以把内部循环用矩阵运算替代但作为作业验证循环版本反而更容易调试和检查中间变量。4.2 抽取、插值与延迟对齐的实操细节抽取和插值这两步看着简单出问题最多也最隐蔽。有几个经验值得记下来。抽取的边界位置。如果你对第1个样本开始抽取会丢掉滤波器的群延迟。最好先对卷积结果做时间对齐再从第0个位置开始抽取。上面的代码里已经先做了(N-1)/2的裁剪这就是在把滤波器的中心延迟消掉。不做这一步系统输出的相位特性会乱频谱图上看不出来但重构的波形时序会整体偏移。插值补零后的增益问题。插值是对子带信号每两个样本之间插入D-1个零这不会改变信号能量但综合滤波器输出后信号的幅度会受到综合滤波器增益的影响。如果原型滤波器没有做归一化M路相加后信号幅度会放大M倍。所以设计阶段一定要对原型滤波器做能量归一化$$h(n) \frac{h(n)}{\sum_k |G_k(e^{j\omega})|^2}$$一个简单做法是在设计原型时除以sum(h)这样各子带在通带内单位增益相加后整体输出增益为1。做完归一化后再测一次重构增益理论上应该接近1误差一般来自过渡带的非理想叠加。延迟对齐是最容易被忽视的。完全重建要求输出是输入的延迟副本。用上述代码里的tau_total截取后你可以直接计算重构误差。但如果延迟算错误差会大得离谱而且你会观察到误差曲线的形状像是一个“回声”而不是噪声。测真实延迟有一个很笨但很可靠的土办法构造一个单位冲激 $x(n) \delta(n)$送入整个链路然后找输出绝对值最大的峰值位置那个位置就是真实延迟。把测量的延迟值替换掉公式计算值重构误差会瞬间降下来。这个方法在你改过滤波器长度、改了抽取因子之后尤其好用。4.3 重构误差的三个度量指标验证完全重建不能只看波形“像不像”。要用数值指标说话。我建议至少算三个指标重构误差能量比dB$$E_{\text{rms}} 10 \log_{10} \left( \frac{\sum_n (x(n) - \hat{x}(n))^2}{\sum_n x(n)^2} \right)$$这个值越小越好。好的QAM滤波器组轻松做到-60dB以下。幅度失真曲线$$A(\omega) |T_0(e^{j\omega})|$$用一组不同频率的正弦信号分别通过链路记录输出与输入的幅度比画出幅度响应曲线。理想情况下是平直线起伏越大说明幅度失真越明显。混叠泄漏指标$$A_{\text{alias}}(\omega) \frac{\sum_{l1}^{D-1} |T_l(e^{j\omega})|^2}{|T_0(e^{j\omega})|^2}$$这个指标用来衡量所有混叠项的总能量相对于正常信号传输的比例。可以拿一段在某个子带之外有明显能量的宽带信号送入系统检查输出种是否出现了输入没有的频率成分。写作业报告的时候给一段代码自动计算这三个指标function report check_reconstruction(x, y, fs) % 计算重构质量指标 % 输入/输出信号对齐后长度一致 % 1. 重构误差能量比 e x - y; err_db 10 * log10(sum(e.^2) / (sum(x.^2) eps)); % 2. 频域幅度失真 X fft(x); Y fft(y); mag_dist abs(Y) ./ (abs(X) eps); % 3. 高频混叠分量检测简单近似 % 对高频段做能量比较正常信号无高频残留时此项应为零 hf_energy sum(abs(Y(round(fs*0.75):end)).^2); total_energy sum(abs(Y).^2); alias_ratio hf_energy / (total_energy eps); report struct(err_db, err_db, ... mag_dist, mag_dist, ... alias_ratio, alias_ratio); fprintf(重构误差能量比: %.2f dB\n, err_db); fprintf(混叠分量占比: %.4f%%\n, alias_ratio * 100); end如果你的系统设计正确误差能量比应该在-60dB以下混叠分量占比应该是10的负几次方级别。如果达不到要么滤波器阻带衰减不够要么延迟没对齐要么调制公式写错了。4.4 与余弦调制滤波器组的关联复调制QAM滤波器组虽然结构清晰但在实际数据压缩系统中更常见的是余弦调制滤波器组CMFB。原因很简单复调制产生的是复数子带信号每个样本由实部和虚部组成数据量直接翻倍而余弦调制只产生实信号不会增加数据量对压缩系统更友好。余弦调制滤波器组的分析滤波器为$$h_k(n) 2h(n) \cos\left(\frac{\pi(k0.5)}{M}(n - \tfrac{N-1}{2}) \theta_k\right)$$其中 $\theta_k$ 是相位参数常见取法为 $\theta_k (-1)^k \cdot \pi/4$。它的本质是把复QAM滤波器组的M个子带“折叠”成2M个实数子带再截取其中M个从而做到临界采样下的完全重建。如果你仔细看公式会发现余弦其实是两个复指数的和所以CMFB可以被视为QAM滤波器组的实数化变形。作业如果要求的是“QAM滤波器组”用复调制路线最直接。但如果你的课程后续有“把子带信号编码成比特流”的环节我建议把CMFB也一并实现因为它输出的实数子带信号可以直接进量化器。实现方式是在复调制的基础上把实部和虚部合并成余弦形式其余流程不变。报告里加上两种结构的对比导师会觉得你有深度。5. 作业中常见的坑排查与解决5.1 重构波形“看起来”对但误差很大我见过不少同学卡在这个坑里跑完系统画出来波形两条曲线肉眼几乎重合但一算误差发现是-20dB级别远达不到完全重建的标准。这种情况通常是“相位对了但幅度不对”或者“幅度对了但混叠没有完全抵消”。排查顺序建议如下首先确认增益。在链路中任意位置插入一个增益缩放不是好事最可能的原因是原型滤波器没有归一化。你可以在综合端输出后用正弦信号测一下输出幅度如果是输入的2倍或者0.5倍那就去检查原型滤波器的能量归一化。其次看延迟。如果误差曲线呈现出差信号像一个“位移后的信号副本”而非噪声几乎可以断定是延迟没对齐。把延迟改成冲激测量值试试。最后看混叠。如果误差形状是高频抖动幅度在几十到几百之间波动那就是混叠泄漏。检查阻带衰减是否够大或者子带截止频率是否设错了。把原型滤波器的频率响应画出来看每个子带之间是否留有保护间隔。没有保护间隔的过采样配置混叠泄漏是很难压下去的。5.2 延迟补偿到底该怎么算这是一个让很多人头大的问题。教科书里的总延迟公式依赖滤波器实现的细节你换了边界处理方式延迟就可能变。最可靠的办法是“冲激定位法”简单粗暴但绝对有效。具体做法输入一个单点冲激 $x [1, 0, 0, \dots]$经过完整的分析-综合链路用find(abs(y) max(abs(y)))找到输出波形的尖峰位置那个位置就是延迟样本数。把这个数代进你的对齐代码里重构误差立刻会变得好看很多。需要注意冲激定位的前提是系统基本正确只是延迟没对齐。如果系统里有实质性的混叠问题冲激响应也会变得一团乱定位到的峰值不一定可靠。所以先用正弦全频带扫一遍确认系统没有致命问题再做延迟对齐排查效率更高。5.3 频域混叠泄漏的定位技巧混叠泄漏定位起来最麻烦的点在于它和想要的信号混在一起时域上根本看不出差别。有一个很好用的方法输入一个窄带正弦信号频率设置为某个子带的通带中心附近另一个子带通带之外。如果系统完全重建成功输出应该只包含这个正弦频率频谱图上只有一个干净的峰。但如果存在混叠泄漏输出频谱会出现一个“鬼峰”位置在其他子带的中心频率上。通过鬼峰出现的位置你可以反推是哪个子带滤波器组没配合好。我调试时常用这段代码% 输入一个窄带信号观察混叠 f_test 0.31; % 归一化频率位于第二个子带附近 n 0:1023; x sin(2 * pi * f_test * n); y analysis_synthesis(x, H, G, D); Xf abs(fft(x)); Yf abs(fft(y)); % 画图比较肉眼找鬼峰 plot(linspace(0, 1, 1024), Xf, LineWidth, 1.5); hold on; plot(linspace(0, 1, 1024), Yf, LineWidth, 1.5); legend(输入,输出);如果输出频带里出现输入没有的频率分量说明混叠没消干净。这时候调整的无非是这几个参数原型滤波器阻带衰减、子带调制相位的偏移量、综合滤波器配对方式。三项逐一排查基本能定位到问题根因。5.4 数值精度与边界处理滤波器组的卷积运算本质上是对浮点数做大量乘加。滤波器越长误差累积越明显。如果用的是单精度浮点重构误差可能被浮点噪声限制在-50dB左右怎么调都过不去。作业环境下建议用双精度MATLAB默认就是双精度一般能跑到-70dB以下满足完全重建的验收。边界处理也很关键。conv函数的默认行为是不截断的卷积结果比原始信号长如果直接截取固定长度头部和尾部会留下因为滤波器非因果部分引起的暂态。要做到真正的完全重建需要把暂态部分丢弃只保留不受边界影响的“稳态区”。前文代码里已经做了这个截取原理是把滤波器中心延迟去掉后再取等长输出。如果你在信号两端做了补零或延拓要注意补零会引入额外的高频瞬态可能导致输出信号首尾几百个样本的误差特别大。解决办法是评测误差时跳过首尾的保护区间只统计中间稳定区间的重构误差。写报告时注明这一点说明你知道边界效应的影响范围属于加分项。6. 最后分享一点个人体会完全重建QAM滤波器组这个作业做完之后我对“信号处理系统的验收”有了完全不同的认识。以前我以为滤波就是卷积、变换就是FFT跑通了就完事。但实际上一个系统是否真正合格要看它在实际链路中的误差表现而不是单个模块的理想假设。这个作业逼着我把每个环节的延迟、增益、相位都算得清清楚楚才算把多速率信号处理的框架真正装进脑子里。如果你也在做这个作业我的建议是别急着写代码先把第2章那三个失真项的公式手推一遍再用代码实现效率会高很多。代码出了问题回到公式去找比盲改参数靠谱得多。调通之后留下那套冲激定位延迟的测试脚本后续做余弦调制滤波器组、MDCT变换器甚至做子带编码都能复用它。最后再分享一个小技巧提交作业时除了重构误差附上一张输入输出信号的频谱叠加图能直观体现完全重建的效果。老师一眼就能看出你的系统有没有混叠泄漏比任何文字说明都管用。
数据压缩课最麻烦的作业往往不是写熵编码或者量化器而是让你把一整套滤波器组从原理到仿真完整实现出来。完全重建QAM滤波器组这个题目表面上是调库调用实际上牵扯到多速率信号处理、原型滤波器设计、调制映射、延迟对齐和混叠分析每个环节都能埋雷。这篇文章把我做这个作业时踩过的坑、总结过的推导、以及可以直接跑通的实现步骤都整理出来给正在被这门作业折磨的人一个参考。1. 作业到底在做什么QAM滤波器组与子带编码1.1 数据压缩作业里的“滤波器组”角色数据压缩的核心思路无非是去掉冗余。时域预测编码去的是时间冗余变换编码去的是空间或统计冗余而滤波器组做的事情是帮你在频域上把信号切成若干个子带然后根据每个子带的能量和感知重要程度分配合适的量化位数。这个过程在语音编码、音频编码和图像压缩里都有应用经典如子带图像编码、MDCT音频编码背后都是滤波器组的思路。你可能会问既然有FFT这种“一刀切”的变换为什么还要用滤波器组区别在于FFT是对整段信号做固定分块块与块之间没有任何平滑过渡会导致很明显的边界效应。而滤波器组的设计天然带有重叠和过渡分析端把信号划分到子带综合端又能把子带重新拼回完整信号中间的切换是连续的重构质量就靠“完全重建”这个性质来保证。作业里提到的QAM滤波器组就是其中一种实现方式。它名字里带着“QAM”容易让人联想到通信里的正交幅度调制实际上它跟通信里的星座图没有直接关系。这里的QAM是指用复指数调制正交调制把原型低通滤波器搬移到不同频带从而形成一组相互覆盖整个频带的子带滤波器。这种设计的好处是结构统一分析端和综合端的滤波器关系明确配合上过采样处理可以做到工程意义上的完全重建。1.2 QAM滤波器组的名词拆解与整体结构先把滤波器的组立结构讲清楚。一个标准的M通道滤波器组包含分析滤波器和综合滤波器两部分。分析端有M个滤波器编号从0到M-1。每个分析滤波器 $H_k(z)$ 的作用是从原始信号中提取第k个子带的成分。原始信号经过这M个滤波器后得到M路信号每一路都只保留了一个频带内的内容然后经过抽取降低采样率。综合端也有M个滤波器编号对应。M路子带信号先经过插值恢复到原始采样率再通过综合滤波器 $G_k(z)$ 滤掉插值产生的镜像谱最后把所有子带输出加起来得到重建信号。QAM滤波器组的结构特点是所有分析滤波器都由同一个低通原型滤波器 $h(n)$ 调制而来。具体做法是$$h_k(n) h(n) \cdot e^{j \omega_k (n - \tfrac{N-1}{2})}$$其中 $\omega_k$ 是第k个子带的中心频率$N$ 是原型滤波器长度。这个“调制”操作本质就是把低通滤波器的频谱搬到 $\omega_k$ 附近。$\omega_k$ 的取法有讲究好的取值能让子带之间刚好无缝衔接也方便后续综合端重建时消除混叠。用这种“一个原型调制”的方式你只需要设计好一个低通滤波器剩下的M个带通滤波器全部由公式生成调试起来很舒服。后续如果要改通道数也不用推翻重来改个M重新生成一遍调制系数就行。这个特性写作业和做工程都很省心。1.3 为什么“完全重建”是硬指标滤波做得好不好最终要落到一个验收标准上原始信号 $x(n)$ 经过分析-综合系统后输出 $\hat{x}(n)$ 与 $x(n)$ 之间差多少。如果系统是理想的那么除了固定延迟和幅度缩放输出应该和输入一模一样。这种“输入输出只差一个延迟和常数增益”的性质就叫完全重建。对数据压缩而言完全重建意味着只要量化环节不做任何有损操作或者量化误差为零整个编解码链路就不会引入额外失真。在这个前提下压缩的误差完全由量化步长决定系统设计师可以放心地调整压缩参数而不需要担心滤波器组本身在偷偷制造噪声。但现实是滤波器组涉及抽取和插值必然会产生混叠。如果原型滤波器设计得不好或者调制参数不对混叠就会泄漏到重构信号中表现出来就是高频噪声或者周期性的波纹。这些问题你在时域里看波形可能不明显但在频域里一看就露馅了。所以判断作业做得对不对先检查重构误差再检查混叠泄漏。这两个指标合格才算真正把完全重建做出来了。2. 完全重建背后的三个失真项2.1 分析-综合系统的输入输出关系要理解完全重建先要把整个链路的数学表达式写出来。这里我们假设抽取因子为 $D$插值因子也是 $D$。输入端信号 $X(z)$ 经过第k个分析滤波器 $H_k(z)$ 后被 $D$ 倍抽取再被 $D$ 倍插值最后经过综合滤波器 $G_k(z)$ 后求和。这个过程中抽取会产生频谱折叠。只要分析滤波器的阻带不是零点抽取后的子带信号中就会混入相邻频带的成分。插值会产生镜像谱需要在综合端通过滤波器抑制。把这些效应都考虑进去后输出可以表示为$$Y(z) \frac{1}{D}\sum_{l0}^{D-1} X(zW_D^l) \sum_{k0}^{M-1} G_k(z) H_k(zW_D^l)$$其中 $W_D e^{-j2\pi/D}$。这个式子看着吓人其实拆开看就两件事求和符号里 $l0$ 的那一项对应的是原始信号的正常传输路径而 $l \neq 0$ 的各项对应的是因为抽取折叠而产生的混叠成分。每一项混叠的强度取决于分析滤波器在折叠频点处的响应和综合滤波器在该处的抑制能力。为了简化表达定义传输函数$$T_l(z) \frac{1}{D}\sum_{k0}^{M-1} G_k(z) H_k(zW_D^l)$$那么输出就是$$Y(z) T_0(z) X(z) \sum_{l1}^{D-1} T_l(z) X(zW_D^l)$$到这里完全重建的三个条件就浮出水面了。2.2 混叠、幅度失真、相位失真分别如何产生混叠失真对应的是 $T_l(z) 0$对所有 $l \neq 0$ 成立。也就是说那些因为抽取而折叠进来的频谱成分必须在综合端被彻底抵消。理想情况下分析滤波器的阻带足够干净抽取后折叠进来的成分本来就很小但更重要的一点是综合滤波器的频响覆盖范围要能补回折叠丢失的部分。在QAM滤波器组中由于各子带滤波器是同一个原型调制得到的相邻子带的尾部响应可以相互抵消这正是调制结构比任意滤波器组的优势所在。幅度失真对应的是 $T_0(z)$ 的幅频响应必须是常数或者近似常数。如果 $T_0(z)$ 在某些频率处有明显的谷或峰那么原信号经过系统后某些频段会被削弱或放大听起来就是“频谱不平坦”。对于完全重建系统要求 $|T_0(e^{j\omega})|$ 在所有频率上保持一致通常归一化为1。相位失真对应的是 $T_0(z)$ 的相位必须与频率成线性关系也就是整体只引入一个固定延迟 $\tau$。如果相位呈现非线性信号的不同频率分量会被错开波形就会产生畸变。线性相位原型滤波器在这里帮了大忙因为只要原型 $h(n)$ 是线性相位的调制后的分析滤波器也保持线性相位特性综合端再取合适配对整体链路的相位就可控。这三个条件都满足时$Y(z) c z^{-\tau} X(z)$也就是完全重建。实际操作中工程意义上的完全重建允许极小的偏差比如重构误差低于-60dB就可以当作“重建成功”来对待。2.3 过采样QAM为何能实现完全重建你可能在教科书里看到过临界采样的M通道滤波器组要做完全重建原型滤波器的设计条件非常苛刻尤其是M比较大的时候几乎很难同时压住混叠和幅度失真。那为什么QAM滤波器组还能做关键在于它可以配合过采样使用。过采样的含义是分析滤波器个数是M但抽取因子D小于M。最常见的配置是 $D M/2$也就是2倍过采样。原本信号被切成M个子带每个子带只保留 $1/M$ 的带宽信息按理说采样率可以降为原来的 $1/M$但过采样只降到 $1/D$等于每个子带信号里留了余量频谱之间留有保护带抽取折叠时混叠就小得多。打个比方临界采样就像把一整条面包切成M片每片都紧挨着切歪一点就露馅。过采样则是让每片之间留一点空隙切成M片但每片都薄一点拼接的时候自然不容易出错。对于完全重建来说保护带就是容错空间。QAM调制再加上过采样综合端重建时那些 $l \neq 0$ 的混叠项会被原型滤波器的阻带衰减和过渡带设计联合压制重构质量可以做得非常高。这也是为什么作业里推荐使用过采样QAM而不是临界采样QMF前者在工程实现上简单得多对原型滤波器的要求也更宽松。理解了这点后面设计参数时你就知道哪些地方可以放松哪些地方必须死磕。3. 原型滤波器设计与调制构造步骤3.1 设计M通道系统的关键参数做QAM滤波器组第一步是定参数。参数之间互相制约定错了后面返工成本很高。我整理了一张参数表写作业和做毕设都可以直接参考。参数符号取值依据我的推荐通道数M子带数量取决于压缩需求M8或M16量级适中抽取因子D过采样率决定 DMDM/2典型2倍过采样原型滤波器长度N越长过渡带越窄但延迟越大N4KM-1K取2~4原型截止频率fc要覆盖子带占用带宽fc1/(2*D)归一化频率阻带衰减AS决定混叠泄漏水平60dB以上最好80dB过渡带宽Δf与N直接相关由AS和N反推这里尤其要注意原型滤波器长度和过渡带宽的关系。窗函数法设计滤波器时过渡带宽与滤波器长度近似成反比。N太小过渡带宽就宽相邻子带会互相污染N太大整个系统的延迟也大而且接近长滤波器的数值精度问题会浮现出来。我的经验是先定M再定D然后用 $N 4KM - 1$ 作为初始值如果重构误差不够再把K往上调。截止频率的取值也别拍脑袋。分析滤波器的通带要覆盖 $\pi/M$ 范围内的信息而抽取因子D决定了子带信号保留的带宽。当 $D M/2$ 时每个子带的保留带宽约为 $1/(2D) 1/M$所以原型滤波器的通带截止频率取 $\pi/M$ 附近比较合理。你可以把原型滤波器的6dB衰减点设置在这个位置保证相邻子带的交叠处能量互补。3.2 用窗函数法设计原型低通滤波器确定好参数后设计原型滤波器最省事的方法就是窗函数法尤其是Kaiser窗因为它允许你独立控制阻带衰减和过渡带宽不像矩形窗那样参数是死的。Kaiser窗设计低通滤波器的步骤如下确定阻带衰减 $A_s$单位dB比如取80dB。根据衰减值计算Kaiser窗的 $\beta$ 参数$$\beta \begin{cases} 0.1102(A_s - 8.7), A_s 50 \ 0.5842(A_s - 21)^{0.4} 0.07886(A_s - 21), 21 \le A_s \le 50 \end{cases}$$确定过渡带宽 $\Delta\omega$然后反推滤波器长度$$N \frac{A_s - 7.95}{2.285 \cdot \Delta\omega} 1$$过渡带宽要根据子带间距来定。子带之间能容忍的过渡带范围大约是 $0.1\pi/M$ 到 $0.4\pi/M$再宽就会侵占相邻子带的频带。用MATLAB做这一步非常直接。下面是我作业里实际用过的设计代码function h design_prototype(M, D, K, As) % 设计Kaiser窗原型低通滤波器 % M: 通道数 % D: 抽取因子 % K: 长度系数滤波器长度为 4*K*M - 1 % As: 阻带衰减 (dB) N 4 * K * M - 1; fc 1 / (2 * D); % 由阻带衰减计算Kaiser窗beta参数 if As 50 beta 0.1102 * (As - 8.7); else beta 0.5842 * (As - 21)^0.4 0.07886 * (As - 21); end % 理想低通脉冲响应线性相位 n 0:N-1; h_ideal 2 * fc * sinc(2 * fc * (n - (N-1)/2)); % 加Kaiser窗 w kaiser(N, beta); h h_ideal .* w; % 归一化使各子带合成增益为常数 h h / sum(h); end这里有一个细节很多教程不提为什么长度取 $4KM-1$ 而不是随便一个偶数因为QAM滤波器组在做调制时需要原型滤波器关于中心点对称即满足 $h(n) h(N-1-n)$这要求N是奇数。同时长度要是M的倍数关系这样调制后的相位关系才规整综合端重建时各子带之间的相位对齐才有保证。$4KM-1$ 这个形式保证了N既是奇数又与M有明确的模关系后续分析延迟和做多相分解都会省事很多。3.3 从原型到QAM滤波器组调制这一步怎么算原型滤波器设计好后剩下的就是调制。QAM滤波器组使用复指数调制第k个子带的中心频率为$$\omega_k \frac{2\pi(k 0.5)}{M}, \quad k 0, 1, \dots, M-1$$这里加了一个 $0.5$ 的偏移不是从0开始而是从半整数倍开始。这样做的目的是让第一个子带的中心频率离开0频避免DC分量对第一个子带产生干扰同时也让原型滤波器在0频附近的响应不必完美因为第一个子带不再覆盖直流点。分析滤波器的生成公式$$h_k(n) h(n) \cdot e^{j\omega_k (n - \tfrac{N-1}{2})}$$注意调制指数里的 $n - (N-1)/2$ 这个平移项。加上它是为了保证 $h_k(n)$ 仍然保持关于中心点的共轭对称性也就是$$h_k(N-1-n) h_k^*(n)$$这个性质在后面做延迟补偿时很有用。写代码的时候平移项千万别漏否则滤波器的相位会变得很乱重构出来的信号波形会歪得离谱。生成所有分析滤波器的MATLAB代码function [H, G] design_qam_filterbank(h, M) % 基于原型滤波器h生成M通道QAM分析滤波器组H和综合滤波器组G N length(h); n (0:N-1) - (N-1)/2; % 中心对齐的时间索引 H zeros(M, N); G zeros(M, N); for k 0:M-1 w_k 2 * pi * (k 0.5) / M; mod exp(1j * w_k * n); H(k1, :) h .* mod; % 综合滤波器取共轭对称形式保证相位对齐 G(k1, :) conj(fliplr(H(k1, :))); end end综合滤波器取分析滤波器的共轭时间反转这个动作后面会解释。这里先记住分析端的滤波器负责把信号拆开综合端的滤波器负责把拆开的信号还原两者是配对关系不是随意组合。3.4 分析/综合滤波器配对与延迟补偿综合滤波器取共轭时间反转是有讲究的。因为原型滤波器是线性相位的调制后的分析滤波器具备共轭对称性。综合滤波器取 $G_k(z) z^{-(N-1)} H_k^(1/z^)$在实数系下等价于对时间索引反转并取共轭。这样做可以保证整个分析-综合链路的传递函数中非延迟项相互抵消只留下一个整体延迟。如果综合滤波器随便取比如直接用分析滤波器本身而不做任何反转那么链路中每个子带的相位变化会叠加输出信号会发生频率相关的相位扭曲。在作业的验证阶段你可能发现重构波形幅度差不多但过零点位置全乱了那就是配对方式出了问题。延迟补偿也是关键。滤波器组整体链路是有通带延迟的对于长度N的原型滤波器分析滤波过程引入的延迟大约是 $(N-1)/2$综合过程同样引入 $(N-1)/2$再加上抽取插值处理带来的处理延迟 $D-1$ 等总的输出延迟大概为$$\tau \frac{N-1}{2} \frac{N-1}{2} (D-1) N - 1 D - 1$$这个公式在不同实现里会有细微差别和你滤波器的边界处理方式有关。最稳妥的办法不是手工推算而是用单位冲激做一次“冲激响应测试”直接测量输出峰值位置得到准确延迟。具体做法在第4章讲。4. 用仿真验证完全重建的完整流程4.1 搭建分析滤波器组与综合滤波器组参数设计完了接下来就是把整个系统跑起来。我建议用一个测试脚本串联整个仿真生成输入信号、通过分析端、过采样处理、通过综合端、对比重构误差。这样改任何参数跑一遍就能看到结果。分析端的处理流程是输入信号 $x(n)$ 分别与M个分析滤波器做卷积得到M路滤波结果然后每路按抽取因子D进行抽取。注意这里抽取要从正确的样本位置开始。MATLAB的downsample函数默认从第1个样本开始但滤波有延迟如果你从第1个样本就抽取等于把最前面的暂态也算进去了重构会不准。正确做法是先滤波再把前 $(N-1)/2$ 个样本丢弃或做补偿然后再抽取。综合端的流程是每路子带信号先按D倍插值补零然后与对应的综合滤波器做卷积最后把M路结果相加。function y analysis_synthesis(x, H, G, D) % 完整的分析-综合链路 % x: 输入信号行向量 % H: 分析滤波器组M x N % G: 综合滤波器组M x N % D: 抽取/插值因子 [M, N] size(H); % 计算总延迟后面用于对齐 tau_total N - 1 D - 1; % 分析端 subbands zeros(M, length(x)); for k 1:M tmp conv(x, H(k, :)); tmp tmp((N-1)/2 1 : end - (N-1)/2); % 去掉滤波暂态 subbands(k, :) tmp; end % 抽取 decimated subbands(:, 1:D:end); % 插值 upsampled zeros(M, size(decimated, 2) * D); for k 1:M upsampled(k, 1:D:end) decimated(k, :); end % 综合端 y zeros(1, size(upsampled, 2) N - 1); for k 1:M tmp conv(upsampled(k, :), G(k, :)); y y tmp; end % 截取有效部分并与输入对齐 y y(tau_total 1 : tau_total length(x)); end这个函数把抽取和插值都显式写出来了方便你观察每一步的形状。实际跑起来如果信号太长可以把内部循环用矩阵运算替代但作为作业验证循环版本反而更容易调试和检查中间变量。4.2 抽取、插值与延迟对齐的实操细节抽取和插值这两步看着简单出问题最多也最隐蔽。有几个经验值得记下来。抽取的边界位置。如果你对第1个样本开始抽取会丢掉滤波器的群延迟。最好先对卷积结果做时间对齐再从第0个位置开始抽取。上面的代码里已经先做了(N-1)/2的裁剪这就是在把滤波器的中心延迟消掉。不做这一步系统输出的相位特性会乱频谱图上看不出来但重构的波形时序会整体偏移。插值补零后的增益问题。插值是对子带信号每两个样本之间插入D-1个零这不会改变信号能量但综合滤波器输出后信号的幅度会受到综合滤波器增益的影响。如果原型滤波器没有做归一化M路相加后信号幅度会放大M倍。所以设计阶段一定要对原型滤波器做能量归一化$$h(n) \frac{h(n)}{\sum_k |G_k(e^{j\omega})|^2}$$一个简单做法是在设计原型时除以sum(h)这样各子带在通带内单位增益相加后整体输出增益为1。做完归一化后再测一次重构增益理论上应该接近1误差一般来自过渡带的非理想叠加。延迟对齐是最容易被忽视的。完全重建要求输出是输入的延迟副本。用上述代码里的tau_total截取后你可以直接计算重构误差。但如果延迟算错误差会大得离谱而且你会观察到误差曲线的形状像是一个“回声”而不是噪声。测真实延迟有一个很笨但很可靠的土办法构造一个单位冲激 $x(n) \delta(n)$送入整个链路然后找输出绝对值最大的峰值位置那个位置就是真实延迟。把测量的延迟值替换掉公式计算值重构误差会瞬间降下来。这个方法在你改过滤波器长度、改了抽取因子之后尤其好用。4.3 重构误差的三个度量指标验证完全重建不能只看波形“像不像”。要用数值指标说话。我建议至少算三个指标重构误差能量比dB$$E_{\text{rms}} 10 \log_{10} \left( \frac{\sum_n (x(n) - \hat{x}(n))^2}{\sum_n x(n)^2} \right)$$这个值越小越好。好的QAM滤波器组轻松做到-60dB以下。幅度失真曲线$$A(\omega) |T_0(e^{j\omega})|$$用一组不同频率的正弦信号分别通过链路记录输出与输入的幅度比画出幅度响应曲线。理想情况下是平直线起伏越大说明幅度失真越明显。混叠泄漏指标$$A_{\text{alias}}(\omega) \frac{\sum_{l1}^{D-1} |T_l(e^{j\omega})|^2}{|T_0(e^{j\omega})|^2}$$这个指标用来衡量所有混叠项的总能量相对于正常信号传输的比例。可以拿一段在某个子带之外有明显能量的宽带信号送入系统检查输出种是否出现了输入没有的频率成分。写作业报告的时候给一段代码自动计算这三个指标function report check_reconstruction(x, y, fs) % 计算重构质量指标 % 输入/输出信号对齐后长度一致 % 1. 重构误差能量比 e x - y; err_db 10 * log10(sum(e.^2) / (sum(x.^2) eps)); % 2. 频域幅度失真 X fft(x); Y fft(y); mag_dist abs(Y) ./ (abs(X) eps); % 3. 高频混叠分量检测简单近似 % 对高频段做能量比较正常信号无高频残留时此项应为零 hf_energy sum(abs(Y(round(fs*0.75):end)).^2); total_energy sum(abs(Y).^2); alias_ratio hf_energy / (total_energy eps); report struct(err_db, err_db, ... mag_dist, mag_dist, ... alias_ratio, alias_ratio); fprintf(重构误差能量比: %.2f dB\n, err_db); fprintf(混叠分量占比: %.4f%%\n, alias_ratio * 100); end如果你的系统设计正确误差能量比应该在-60dB以下混叠分量占比应该是10的负几次方级别。如果达不到要么滤波器阻带衰减不够要么延迟没对齐要么调制公式写错了。4.4 与余弦调制滤波器组的关联复调制QAM滤波器组虽然结构清晰但在实际数据压缩系统中更常见的是余弦调制滤波器组CMFB。原因很简单复调制产生的是复数子带信号每个样本由实部和虚部组成数据量直接翻倍而余弦调制只产生实信号不会增加数据量对压缩系统更友好。余弦调制滤波器组的分析滤波器为$$h_k(n) 2h(n) \cos\left(\frac{\pi(k0.5)}{M}(n - \tfrac{N-1}{2}) \theta_k\right)$$其中 $\theta_k$ 是相位参数常见取法为 $\theta_k (-1)^k \cdot \pi/4$。它的本质是把复QAM滤波器组的M个子带“折叠”成2M个实数子带再截取其中M个从而做到临界采样下的完全重建。如果你仔细看公式会发现余弦其实是两个复指数的和所以CMFB可以被视为QAM滤波器组的实数化变形。作业如果要求的是“QAM滤波器组”用复调制路线最直接。但如果你的课程后续有“把子带信号编码成比特流”的环节我建议把CMFB也一并实现因为它输出的实数子带信号可以直接进量化器。实现方式是在复调制的基础上把实部和虚部合并成余弦形式其余流程不变。报告里加上两种结构的对比导师会觉得你有深度。5. 作业中常见的坑排查与解决5.1 重构波形“看起来”对但误差很大我见过不少同学卡在这个坑里跑完系统画出来波形两条曲线肉眼几乎重合但一算误差发现是-20dB级别远达不到完全重建的标准。这种情况通常是“相位对了但幅度不对”或者“幅度对了但混叠没有完全抵消”。排查顺序建议如下首先确认增益。在链路中任意位置插入一个增益缩放不是好事最可能的原因是原型滤波器没有归一化。你可以在综合端输出后用正弦信号测一下输出幅度如果是输入的2倍或者0.5倍那就去检查原型滤波器的能量归一化。其次看延迟。如果误差曲线呈现出差信号像一个“位移后的信号副本”而非噪声几乎可以断定是延迟没对齐。把延迟改成冲激测量值试试。最后看混叠。如果误差形状是高频抖动幅度在几十到几百之间波动那就是混叠泄漏。检查阻带衰减是否够大或者子带截止频率是否设错了。把原型滤波器的频率响应画出来看每个子带之间是否留有保护间隔。没有保护间隔的过采样配置混叠泄漏是很难压下去的。5.2 延迟补偿到底该怎么算这是一个让很多人头大的问题。教科书里的总延迟公式依赖滤波器实现的细节你换了边界处理方式延迟就可能变。最可靠的办法是“冲激定位法”简单粗暴但绝对有效。具体做法输入一个单点冲激 $x [1, 0, 0, \dots]$经过完整的分析-综合链路用find(abs(y) max(abs(y)))找到输出波形的尖峰位置那个位置就是延迟样本数。把这个数代进你的对齐代码里重构误差立刻会变得好看很多。需要注意冲激定位的前提是系统基本正确只是延迟没对齐。如果系统里有实质性的混叠问题冲激响应也会变得一团乱定位到的峰值不一定可靠。所以先用正弦全频带扫一遍确认系统没有致命问题再做延迟对齐排查效率更高。5.3 频域混叠泄漏的定位技巧混叠泄漏定位起来最麻烦的点在于它和想要的信号混在一起时域上根本看不出差别。有一个很好用的方法输入一个窄带正弦信号频率设置为某个子带的通带中心附近另一个子带通带之外。如果系统完全重建成功输出应该只包含这个正弦频率频谱图上只有一个干净的峰。但如果存在混叠泄漏输出频谱会出现一个“鬼峰”位置在其他子带的中心频率上。通过鬼峰出现的位置你可以反推是哪个子带滤波器组没配合好。我调试时常用这段代码% 输入一个窄带信号观察混叠 f_test 0.31; % 归一化频率位于第二个子带附近 n 0:1023; x sin(2 * pi * f_test * n); y analysis_synthesis(x, H, G, D); Xf abs(fft(x)); Yf abs(fft(y)); % 画图比较肉眼找鬼峰 plot(linspace(0, 1, 1024), Xf, LineWidth, 1.5); hold on; plot(linspace(0, 1, 1024), Yf, LineWidth, 1.5); legend(输入,输出);如果输出频带里出现输入没有的频率分量说明混叠没消干净。这时候调整的无非是这几个参数原型滤波器阻带衰减、子带调制相位的偏移量、综合滤波器配对方式。三项逐一排查基本能定位到问题根因。5.4 数值精度与边界处理滤波器组的卷积运算本质上是对浮点数做大量乘加。滤波器越长误差累积越明显。如果用的是单精度浮点重构误差可能被浮点噪声限制在-50dB左右怎么调都过不去。作业环境下建议用双精度MATLAB默认就是双精度一般能跑到-70dB以下满足完全重建的验收。边界处理也很关键。conv函数的默认行为是不截断的卷积结果比原始信号长如果直接截取固定长度头部和尾部会留下因为滤波器非因果部分引起的暂态。要做到真正的完全重建需要把暂态部分丢弃只保留不受边界影响的“稳态区”。前文代码里已经做了这个截取原理是把滤波器中心延迟去掉后再取等长输出。如果你在信号两端做了补零或延拓要注意补零会引入额外的高频瞬态可能导致输出信号首尾几百个样本的误差特别大。解决办法是评测误差时跳过首尾的保护区间只统计中间稳定区间的重构误差。写报告时注明这一点说明你知道边界效应的影响范围属于加分项。6. 最后分享一点个人体会完全重建QAM滤波器组这个作业做完之后我对“信号处理系统的验收”有了完全不同的认识。以前我以为滤波就是卷积、变换就是FFT跑通了就完事。但实际上一个系统是否真正合格要看它在实际链路中的误差表现而不是单个模块的理想假设。这个作业逼着我把每个环节的延迟、增益、相位都算得清清楚楚才算把多速率信号处理的框架真正装进脑子里。如果你也在做这个作业我的建议是别急着写代码先把第2章那三个失真项的公式手推一遍再用代码实现效率会高很多。代码出了问题回到公式去找比盲改参数靠谱得多。调通之后留下那套冲激定位延迟的测试脚本后续做余弦调制滤波器组、MDCT变换器甚至做子带编码都能复用它。最后再分享一个小技巧提交作业时除了重构误差附上一张输入输出信号的频谱叠加图能直观体现完全重建的效果。老师一眼就能看出你的系统有没有混叠泄漏比任何文字说明都管用。
