多速率信号处理与采样率转换:抗混叠滤波器设计实战指南
2026/9/15 4:28:35 网站建设 项目流程

简介:面向数字信号处理学习者与MATLAB开发者,这份4KB的小型代码包聚焦多速率信号处理中的采样率转换与抗混叠滤波器设计。资源共2个文件,包含可运行的main.m脚本与配套README.md说明文档。脚本基于MATLAB信号处理工具箱完成信号插值、抽取及抗混叠低通滤波的仿真,配合文档可直观理解resample、upfirdn等核心函数的触发条件、参数设置与结果验证方法,并掌握不同采样率转换场景下滤波器阶数与截止频率的选取思路。包体虽小,但结构紧凑,适合作为课程设计、毕业设计或工程入门实践的参考模板,帮助读者在短时间内搭建完整的采样率转换仿真流程,并延伸到数字通信、语音处理等实际应用。目前已有65人学习下载,可作为信号处理方向快速上手的轻量示例。

1. 多速率信号处理与采样率转换:为什么设计系统先设计抗混叠滤波器

通信基站的ADC输出采样率是65.536 MHz,符号速率却是2 Msps;音频设备要在48 kHz与44.1 kHz之间来回转换;雷达数字中频之后还要做十倍以上的抽取。这些场景指向同一个问题:采样率转换。多数人的第一反应是抽点或插值,但工程上做完之后的结论通常是反过来的——先算清楚频谱折叠发生在哪个频点,再决定抗混叠滤波器长什么样,否则所有指标都是空谈。

多速率信号处理的核心不是“多采样率”本身,而是速率变化前后频谱如何搬移、叠加、镜像。抽取会把高频折叠进低频,插值会在通带外侧制造镜像副本,这两种伪影都必须靠预滤波器压到目标底噪以下。本文从抽取与插值的频率规则讲起,用MATLAB把几种常用采样率转换路径跑通,再落到抗混叠滤波器的参数设计与验证上。内容适合刚接触数字信号处理的工程师,也适合那些常年用resample但说不清边界条件的老手。

2. 采样率转换的频谱折叠原理:抽取、插值在MATLAB里如何验证

2.1 抽取的频谱扩张:折叠进带内的信号来自哪一段

抽取就是每M点保留一点,x_down[n] = x[Mn]。频域对应关系为X_d(e^{jω}) = 1/M Σ X(e^{j(ω−2πk)/M}),含义是原频谱被拉伸到M倍,然后以2π为周期做叠加。当原始信号中最高频率成分超过fs/(2M)时,超出部分就会折叠回低频区,和有效信号叠在一起。这个折叠不是噪声,是确定性的频谱搬移,所以混叠具体出现在哪个频点是可计算的。

用MATLAB验证这一段规则最直接。生成1 kHz与2.5 kHz两个正弦,16 kHz采样,M=4直接抽取。抽取后采样率变成4 kHz,奈奎斯特边界只有2 kHz,2.5 kHz会折叠到4 kHz − 2.5 kHz = 1.5 kHz,而1 kHz的位置不变。

fs0 = 16000; t = 0:1/fs0:0.2; x = sin(2*pi*1000*t) + sin(2*pi*2500*t); M = 4; xd = x(1:M:end); % 直接抽取,没有抗混叠滤波 Nfft = 2048; f0 = (0:Nfft/2-1)/Nfft*fs0; fd = (0:Nfft/2-1)/Nfft*fs0/M; plot(f0, 20*log10(abs(fft(x, Nfft)))); hold on; plot(fd, 20*log10(abs(fft(xd, Nfft))), 'r'); legend('原信号 16 kHz', '直接抽取 4 kHz'); xlabel('Hz'); ylabel('dB');

这段代码里fft长度取了2048,直接用正弦做频谱会有泄漏,但峰值位置足够说明问题。运行后能看到原始2.5 kHz的谱峰在抽取后移动到1.5 kHz,而1 kHz谱峰不变,这就是没有抗混叠滤波时混叠发生的典型形态。反过来想:如果有用信号恰好落在1.5 kHz,那么来自2.5 kHz的混叠分量会直接叠加在它上面,后续任何处理都无法再分开。

2.2 插值形成的镜像频谱:低通截止频率怎么定

插值是在每两个原始样本之间插入L−1个零,数学形式是x_up[n] = x[n/L](当n是L的倍数时),频域表现为频谱压缩到原来的1/L,并且在原采样率fs的整数倍位置上复制出镜像副本。以16 kHz采样、1 kHz信号为例,L=4插值后采样率变为64 kHz,镜像出现在1 kHz、15 kHz、17 kHz、31 kHz等位置,其中15 kHz是−1 kHz的副本。

L = 4; xu = upsample(x, L); % 中间插零,不滤波 f_up = (0:Nfft/2-1)/Nfft*fs0*L; plot(f_up, 20*log10(abs(fft(xu, Nfft)))); xlabel('Hz'); ylabel('dB');

插值后的目标是把原始带宽之外的镜像滤掉。低通截止频率取fs/(2L)=2 kHz,就能保留1 kHz成分,并把15 kHz、17 kHz等副本全部压掉。这里的重点是:插值本身不会损坏信号,损坏发生在你忘记滤波的时候。真正决定重建质量的是滤波器通带平坦度和对第一镜像的抑制量。

2.3 抗混叠滤波器的位置:先滤波后抽取,先插值后滤波

把两节结论放在一起:抽取之前必须低通,把超过目标奈奎斯特边界的分量提前去掉;插值之后必须低通,把零插入产生的镜像清掉。滤波器放在这两处不是约定俗成,而是频谱规则决定的。后面章节里出现的resample、upfirdn、firdecim等函数,本质上都是在做“滤波+抽取”或“插值+滤波”的组合。理解了这个顺序,才能看懂MATLAB函数内部为什么要先设计一个多相抗混叠滤波器。

提示:抽取的滤波和插值的滤波虽然都叫抗混叠,但一个针对折叠,一个针对镜像。两者对阻带起始频率的要求不同,设计滤波器时不要套用同一组参数。

3. 用MATLAB实现采样率转换:resample、upfirdn与整数倍抽取/插值函数

3.1 最快路径:用resample完成有理数采样率转换

MATLAB里最常被调用的多速率函数是resample,一句话写法是y = resample(x, p, q),把输入按p/q倍的速率重采样。p和q是互质整数,目标采样率与源采样率之比等于p/q。例如48 kHz转44.1 kHz,比值是0.91875,用rat函数得到分子分母:

fs_in = 48000; fs_out = 44100; [p, q] = rat(fs_out / fs_in); % p=147, q=160 y = resample(x, p, q);

resample内部采用Kaiser窗设计一个多相抗混叠滤波器,然后依次做p倍插值、滤波、q倍抽取。对大多数原型验证场景,这一条命令就够用,不需要自己设计滤波器。resample有两个可选参数值得注意:L控制滤波器长度,L越大过渡带越窄;beta控制Kaiser窗形状,间接决定阻带衰减。默认值在多数场合表现尚可,但如果你对混叠底噪有明确指标,就应该自己设计h并改用3.2节的方式。

注意:resample输入输出长度会做自动对齐,输出长度大致为ceil(length(x)*p/q),但首尾会有瞬态,做符号同步或时延敏感处理时不能忽略。

3.2 自定义滤波器:用upfirdn把插值、滤波、抽取拆开

当你需要精确控制抗混叠滤波器指标时,resample的黑盒就不够用了。upfirdn(x, h, p, q)是更底层的原语:先对x做p倍零插值,与h卷积,再做q倍抽取。三个操作在一条调用里完成,但h需要你预先设计好。resample内部做的事情和它基本一致。

p = 147; q = 160; fs_mid = fs_in * p; % 插值后的中间采样率 7.056 MHz fp = 20000; % 通带边界 fst = 22500; % 阻带起始,需覆盖输出Nyquist df = [0 fp fst fs_mid/2] / (fs_mid/2); % 以中间采样率Nyquist做归一化 a = [1 1 0 0]; % 通带1,阻带0 dev = [1e-4 1e-5]; % 通带纹波与阻带衰减(线性幅度) [n, fo, ao, w] = firpmord(df, a, dev); h = firpm(n, fo, ao, w); y = upfirdn(x, h, p, q);

这里最容易被忽略的是h的频率轴。upfirdn中h工作在插值后的采样率fs_mid上,滤波器设计必须换算到fs_mid的奈奎斯特频率。如果直接拿输入采样率的Nyquist做归一化,实际滤波器通带会宽p倍,镜像一个都滤不掉。这段代码里的fst取22500 Hz,略高于输出奈奎斯特22050 Hz,给过渡带留了少量余量,工程上更稳妥。

3.3 整数倍抽取与插值:firdecim、firinterp和intfilt

当转换比为整数时,firdecim和firinterp是两个专用入口。firdecim(h, M, x)等价于先filter再downsample,但内部对h做多相分解,能显著减少计算量。firinterp同理,做L倍插值后滤波。它们的典型调用如下:

M = 4; h = fir1(256, 1/M); % 截止在fs/(2M)的简单低通 yd = firdecim(h, M, x); % 抽取,输出采样率fs/M L = 4; h2 = fir1(256, 1/L); yu = firinterp(h2, L, x); % 插值,输出采样率fs*L

intfilt是更老的函数,用于生成整数倍插值滤波器,现在大部分新代码已经被firinterp或designMultirateFIR取代。新版本MATLAB里designMultirateFIR更值得关注,它把抗混叠滤波器设计与抽取/插值参数统一成一个接口,R2022b之后我一般在正式项目里优先用它。

3.4 函数选型对比:这四种方式别再混着用

函数完成操作关键参数典型场景
resample插值+滤波+抽取,有理数转换p, q, L, beta快速原型、算法验证
upfirdn插值+滤波+抽取,可自定义hh, p, q需要指定滤波器指标
firdecim滤波+M倍抽取h, M固定整数倍抽取
firinterpL倍插值+滤波h, L固定整数倍插值
designMultirateFIR设计抽取/插值滤波器filter, decim, dev, w新工程推荐

很多人犯的错是用downsample直接做抽取:x(1:M:end)得到的结果没有经过任何滤波,高频混叠全部落在带内。对比一下误用和正确路径的底噪,结论非常明显:

y_err = x(1:4:end); % 直接抽取,无滤波 y_ok = resample(x, 1, 4); % 完整抗混叠链路 snr_err = snr(y_err); % 通常只有20~30 dB snr_ok = snr(y_ok); % 能到60 dB以上

snr函数在Signal Processing Toolbox里,只做指标对比可以接受。真正做音频或通信系统时,还要查看混叠落入带内的具体频点,而不仅仅是总SNR。

4. 抗混叠滤波器设计:窗函数、等波纹、多级级联与48 kHz到44.1 kHz实例

4.1 滤波器指标怎么从系统需求推导出来

设计抗混叠滤波器之前先定三个数:通带边界fp、阻带起始频率fst、阻带衰减As。fp通常是信号本身的有效带宽,fst取决于目标采样率下的奈奎斯特边界或镜像/混叠区最低频率,As由系统底噪决定。以48 kHz到44.1 kHz音频转换为例,信号内容只保留到20 kHz,输出奈奎斯特是22.05 kHz,所以fp=20 kHz、fst=22.05 kHz,阻带衰减通常要求90 dB以上,否则可听底噪里会混入镜像。

设计方法优点缺点适用场景
窗函数法 fir1简单、稳定、线性相位通带纹波与阻带衰减互相牵扯快速滤波、指标不苛刻
等波纹 firpm相同阶数阻带抑制更均衡阶数估计不准时需要迭代抽取/插值混叠抑制
最小二乘 firls通带纹波总体最低阻带可能出现个别超差对总失真敏感的音频链路

4.2 等波纹设计:firpmord配合firpm的完整流程

等波纹法的优势是阻带内抑制均匀,不会出现窗函数法那种靠近截止频率处凹陷更深、远处回升的情况。对混叠抑制来说,所需的是整个阻带都低于某个阈值,等波纹法更贴合这个需求。用firpmord先估阶数,再交给firpm设计,是标准流程。

% 中间采样率 7.056 MHz,归一化以插值后Nyquist为1 fs_mid = 48000 * 147; % 7.056 MHz f_norm = [0 20000 22050 fs_mid/2] / (fs_mid/2); a = [1 1 0 0]; dev = [0.0005 1e-5]; % 通带纹波约0.0087 dB,阻带100 dB [n, fo, ao, w] = firpmord(f_norm, a, dev); h = firpm(n, fo, ao, w);

这段代码会算出非常高的阶数,因为过渡带只有2050 Hz,而中间采样率高达7.056 MHz,相对过渡带宽度不到0.15%。这是多速率设计里的经典两难:转换比越大,镜像频点越密集,对滤波器的要求越高。实际工程中不会让过渡带卡这么紧,一般会把fst放宽到24 kHz甚至更高,让滤波器阶数从数千降到几百,代价是损失一点高频带宽。

4.3 多级级联:什么时候拆成两级更划算

当抽取因子D超过10,并且过渡带相对窄时,单级滤波器阶数会爆炸。把D拆成D1×D2,先按D1做第一次抽取,再用D2做第二次,每一级的过渡带都可以放宽,总滤波器长度通常只有单级的几分之一。经验法则:每级抽取因子尽量控制在10以内,过渡带宽度至少占该级奈奎斯特带宽的20%到50%,否则级联的优势不明显。

% 单级抽取 D=12 h12 = designMultirateFIR(1, 12, 64, 90); y1 = upfirdn(x, h12, 1, 12); % 两级抽取 3×4 h3 = designMultirateFIR(1, 3, 64, 90); h4 = designMultirateFIR(1, 4, 64, 90); y2 = upfirdn(upfirdn(x, h3, 1, 3), h4, 1, 4);

4.4 完整示例:48 kHz音频转44.1 kHz的自定义滤波器链路

把前面的内容串成一个可运行脚本。信号源用20 kHz以内的多音合成,方便在频谱上直接看镜像是否被压住:

fs_in = 48000; t = 0:1/fs_in:1; x = sin(2*pi*1000*t) + 0.5*sin(2*pi*10000*t) + 0.2*sin(2*pi*19900*t); p = 147; q = 160; fs_mid = fs_in * p; f_norm = [0 20000 22050 fs_mid/2] / (fs_mid/2); a = [1 1 0 0]; dev = [0.0005 1e-5]; [n, fo, ao, w] = firpmord(f_norm, a, dev); n = n + rem(n, 2); % 保证线性相位为偶数阶 h = firpm(n, fo, ao, w); y = upfirdn(x, h, p, q); y = y(1:round(length(x)*p/q)); % 截断到预期输出长度 fs_out = fs_in * p / q;

脚本中把滤波器阶数向上取偶,是为了让线性相位群延迟正好是整数样本点,后面做时延补偿时会省事。运行后输出采样率是44100 Hz,频谱图在20 kHz以上应该看不到超过−90 dB的镜像残留。如果看到明显的20 kHz以上谱线,先查f_norm是不是按fs_mid归一化,这是最容易出错的点。

5. 多速率系统的多相分解与验证:群延迟补偿和混叠残留定位

5.1 把滤波器拆开:多相分解到底省在哪里

抗混叠滤波器通常动辄上百阶,直接在整个采样率上做卷积,计算量很高。多相分解的思路是把长度N的FIR滤波器按抽取因子M拆成M个子滤波器,每个子滤波器只在低速率的某一路并行工作。数学形式是H(z) = Σ z^−k E_k(z^M),其中E_k(z^M)对应第k个子多相分量。MATLAB的firdecim和resample内部都是这种结构。

手动验证多相结构的等价性可以用一段很短的MATLAB代码:

M = 4; h = fir1(96, 1/M); hp = reshape([h(:); zeros(M - mod(length(h), M), 1)], length(h)+M-mod(length(h),M), []); % hp的每一列就是一个多相子滤波器

hp的每一列就是一路低速子滤波器,输入序列按相位分配到各路,卷积后再合并。子滤波器长度大致是主滤波器的1/M,每路乘法次数大幅减少。这个结构在多速率系统里几乎是标配,理解它才能理解为什么resample处理长序列时比“先filter后downsample”快那么多。

5.2 群延迟怎么补:线性相位滤波器的延迟是N/2

线性相位FIR滤波器的群延迟是(N−1)/2个样本(以滤波器工作采样率计)。使用upfirdn时,这个延迟折算到输出采样率要除以q。不补偿的话,重采样后的信号整体偏移几十个样本,做符号同步或者数据对齐时会直接错位。

delay_mid = (length(h) - 1) / 2; % 中间采样率下的延迟 delay_out = delay_mid / q; % 换算到输出采样率 y = upfirdn(x, h, p, q); y = y(round(delay_out) + 1 : end); % 去掉瞬态并补偿群延迟

这段代码放在设计流程最后,比用xcorr估计延迟更精确,因为群延迟是已知的,不需要用信号相关性去猜。复数调制信号或最小相位滤波器不适用这条规律,最小相位滤波器的群延迟需要用grpdelay逐点计算。

5.3 用pwelch和往返测试定位混叠残留

验证抗混叠效果不能只看时域波形,把转换前后的功率谱放在同一张图里对比是标准做法。pwelch返回的功率谱密度可以直接看到残留镜像落在哪个频段:

[pxx, f] = pwelch(y, hann(4096), 2048, 8192, fs_out); [pxo, fo] = pwelch(x, hann(4096), 2048, 8192, fs_in); semilogy(fo, sqrt(pxo), 'k'); hold on; semilogy(f, sqrt(pxx), 'r'); xline(22050/44100*2, '--', 'Nyquist');

更严格的做法是往返测试:把信号从48 kHz转到44.1 kHz再转回48 kHz,与原始信号比对。往返转换会同时暴露出混叠、镜像和滤波器非理想过渡带三类问题。若往返后的SNR低于系统指标,把两个方向的滤波器阻带衰减各提高10 dB,再测一次,通常能找到问题出在哪一级。多速率系统的交付验证,多数时候就靠这两步:一张频谱图,一个往返SNR。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询