多相滤波器

抽取与抗混叠 FIR 滤波器

抽取的作用是仅保留每 M 个样本中的 1 个样本，从而降低采样率。

如果要从抽取到，首先需要经过一个低通滤波器，滤除所有高于的信号，从而满足 Nyquist 采样定理。

MATLAB 仿真程序如下（代码由 Claude Sonnet 4.5 生成）

close all; clear; clc;

%% ==================== 系统参数 ====================
fs = 40e6;             % 采样率 40 MHz
dsmp_rate = 4;         % 抽取因子
fs_dec = fs / dsmp_rate;

fprintf('原始采样率 fs = %.2f MHz\n', fs/1e6);
fprintf('抽取后采样率 fs'' = %.2f MHz\n', fs_dec/1e6);
fprintf('抽取后奈奎斯特频率 fs''/2 = %.2f MHz\n', fs_dec/2/1e6);

%% ==================== 多音信号生成 ====================
% 设计4个频率分量，其中2个会在抽取后产生混叠
f1 = 2e6;    % 2 MHz - 在新奈奎斯特范围内
f2 = 4e6;    % 4 MHz - 在新奈奎斯特范围内
f3 = 8e6;    % 8 MHz - 接近新奈奎斯特边界
f4 = 13e6;   % 13 MHz - 超出新奈奎斯特范围，会产生混叠

duration = 0.001;  % 信号时长 1ms
t = 0 : 1/fs : duration - 1/fs;

% 生成多音信号
baseband_tx = sin(2*pi*f1*t) + ...
              0.8*sin(2*pi*f2*t) + ...
              0.6*sin(2*pi*f3*t) + ...
              0.5*sin(2*pi*f4*t);

fprintf('\n信号频率分量:\n');
fprintf('  f1 = %.1f MHz (无混叠)\n', f1/1e6);
fprintf('  f2 = %.1f MHz (无混叠)\n', f2/1e6);
fprintf('  f3 = %.1f MHz (边界)\n', f3/1e6);
fprintf('  f4 = %.1f MHz (会混叠到 %.1f MHz)\n', f4/1e6, abs(f4 - fs_dec)/1e6);

%% ==================== 方案 1: 直接抽取（无 AA 滤波）====================
baseband_tx_dec = baseband_tx(1:dsmp_rate:end);
t_dec = t(1:dsmp_rate:end);

%% ==================== 方案 2: 抗混叠滤波 + 抽取 ====================
% 设计 AA 滤波器：截止频率 = fs_dec/2
fc_aa = fs_dec / 2;
Wn = fc_aa / (fs / 2);  % 归一化截止频率
aa_order = 128;         % 增加阶数以获得更陡峭的滚降

h_aa = fir1(aa_order, Wn);
baseband_tx_aa = filter(h_aa, 1, baseband_tx);
% 补偿群延迟
delay_aa = aa_order / 2;
baseband_tx_aa = baseband_tx_aa(delay_aa+1 : end);
% 对齐时间轴
t_aa = t(delay_aa+1 : end);
% 抽取
baseband_tx_aa_dec = baseband_tx_aa(1:dsmp_rate:end);
t_aa_dec = t_aa(1:dsmp_rate:end);

%% ==================== 频谱分析参数 ====================
N = 65536;  % FFT 点数

%% ==================== 图 1: AA 滤波器频率响应 ====================
figure('Name', 'Anti-Aliasing Filter Design');

subplot(1,2,1);
stem(0:aa_order, h_aa, 'filled', 'MarkerSize', 3);
xlabel('样本索引'); ylabel('幅度');
title(sprintf('AA 滤波器冲激响应 (FIR, order=%d)', aa_order));
grid on;

subplot(1,2,2);
[H_aa, f_aa] = freqz(h_aa, 1, N, fs);
plot(f_aa/1e6, 20*log10(abs(H_aa)), 'b', 'LineWidth', 1.5); hold on;
xline(fc_aa/1e6, 'r--', 'LineWidth', 1.5);
xline([f1 f2 f3 f4]/1e6, 'g:', 'LineWidth', 1);
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title('AA 滤波器频率响应');
legend('滤波器响应', sprintf('截止频率 f_c = %.1f MHz', fc_aa/1e6), ...
       '信号频率分量', 'Location', 'southwest');
grid on; xlim([0 fs/2/1e6]); ylim([-100 5]);
hold off;

%% ==================== 图 2: 四组对比（时域 + 频域）====================
figure('Name', 'Signal Comparison: Aliasing vs Anti-Aliasing');

% --- 子图 1: 原始信号时域 ---
subplot(4, 2, 1);
plot(t(1:800)*1e6, baseband_tx(1:800), 'LineWidth', 1);
xlabel('时间 (μs)'); ylabel('幅度');
title('原始多音信号（时域）');
grid on;

% --- 子图 2: 原始信号频域 ---
subplot(4, 2, 2);
f = (-fs/2 : fs/N : fs/2-fs/N) / 1e6;
X = fftshift(fft(baseband_tx, N)) / length(baseband_tx);
spec_orig = 20 * log10(abs(X) + eps);  % 幅度谱 (dB)
plot(f, spec_orig, 'LineWidth', 1); hold on;
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title(sprintf('原始信号频谱 (f_s = %.1f MHz)', fs/1e6));
grid on; xlim([-20 20]);
yl = ylim;

% 标注各频率分量
xline([f1 f2 f3 f4]/1e6, 'r--', 'LineWidth', 1);
xline(-[f1 f2 f3 f4]/1e6, 'r--', 'LineWidth', 1);
text(f1/1e6, yl(2)-10, sprintf('%.0f', f1/1e6), 'Color', 'r', 'FontSize', 8);
text(f2/1e6, yl(2)-10, sprintf('%.0f', f2/1e6), 'Color', 'r', 'FontSize', 8);
text(f3/1e6, yl(2)-10, sprintf('%.0f', f3/1e6), 'Color', 'r', 'FontSize', 8);
text(f4/1e6, yl(2)-10, sprintf('%.0f', f4/1e6), 'Color', 'r', 'FontSize', 8);
hold off;

% --- 子图 3: 直接抽取时域 ---
subplot(4, 2, 3);
plot(t_dec(1:200)*1e6, baseband_tx_dec(1:200), 'o-', 'LineWidth', 1, 'MarkerSize', 3);
xlabel('时间 (μs)'); ylabel('幅度');
title(sprintf('直接 %d 倍抽取（无 AA 滤波）', dsmp_rate));
grid on;

% --- 子图 4: 直接抽取频域（显示混叠）---
subplot(4, 2, 4);
f_dec = (-fs_dec/2 : fs_dec/N : fs_dec/2-fs_dec/N) / 1e6;
X_dec = fftshift(fft(baseband_tx_dec, N)) / length(baseband_tx_dec);
spec_dec = 20 * log10(abs(X_dec) + eps);
plot(f_dec, spec_dec, 'LineWidth', 1.5); hold on;
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title(sprintf('直接抽取频谱 (f_s'' = %.1f MHz) - 出现混叠', fs_dec/1e6));
grid on; xlim([-6 6]);
yl = ylim;

% 标注新奈奎斯特边界
xline(-fs_dec/2/1e6, 'm--', 'LineWidth', 2);
xline( fs_dec/2/1e6, 'm--', 'LineWidth', 2);
text(fs_dec/2/1e6-0.5, yl(2)-5, sprintf('奈奎斯特\n边界'), ...
    'Color', 'm', 'FontSize', 9, 'FontWeight', 'bold');

% 标注混叠的频率分量
f4_alias = abs(f4 - fs_dec) / 1e6;  % 12MHz 混叠到 2MHz
text(f4_alias, yl(2)-15, sprintf('混叠!\n(%.0f→%.0f MHz)', f4/1e6, f4_alias), ...
    'Color', 'r', 'FontWeight', 'bold', 'FontSize', 9, ...
    'HorizontalAlignment', 'center', 'BackgroundColor', [1 0.9 0.9]);
hold off;

% --- 子图 5: AA 滤波后（抽取前）时域 ---
subplot(4, 2, 5);
plot(t_aa(1:800)*1e6, baseband_tx_aa(1:800), 'LineWidth', 1);
xlabel('时间 (μs)'); ylabel('幅度');
title('AA 滤波后（抽取前）');
grid on;

% --- 子图 6: AA 滤波后（抽取前）频域 ---
subplot(4, 2, 6);
X_aa = fftshift(fft(baseband_tx_aa, N)) / length(baseband_tx_aa);
spec_aa = 20 * log10(abs(X_aa) + eps);
plot(f, spec_aa, 'LineWidth', 1); hold on;
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title('AA 滤波后频谱（抽取前）');
grid on; xlim([-20 20]);
yl = ylim;

xline(-fc_aa/1e6, 'm--', 'LineWidth', 1.5);
xline( fc_aa/1e6, 'm--', 'LineWidth', 1.5);
text(fc_aa/1e6+1, yl(2)-5, sprintf('截止\n%.1f MHz', fc_aa/1e6), ...
    'Color', 'm', 'FontSize', 9, 'FontWeight', 'bold');

% 标注被抑制的高频分量
text(f4/1e6, yl(1)+15, sprintf('%.0f MHz\n已抑制', f4/1e6), ...
    'Color', 'g', 'FontWeight', 'bold', 'FontSize', 9, ...
    'HorizontalAlignment', 'center', 'BackgroundColor', [0.9 1 0.9]);
hold off;

% --- 子图 7: AA 滤波 + 抽取时域 ---
subplot(4, 2, 7);
plot(t_aa_dec(1:200)*1e6, baseband_tx_aa_dec(1:200), 'o-', 'LineWidth', 1, 'MarkerSize', 3);
xlabel('时间 (μs)'); ylabel('幅度');
title(sprintf('AA 滤波 + %d 倍抽取', dsmp_rate));
grid on;

% --- 子图 8: AA 滤波 + 抽取频域（无混叠）---
subplot(4, 2, 8);
X_aa_dec = fftshift(fft(baseband_tx_aa_dec, N)) / length(baseband_tx_aa_dec);
spec_aa_dec = 20 * log10(abs(X_aa_dec) + eps);
plot(f_dec, spec_aa_dec, 'LineWidth', 1.5); hold on;
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title(sprintf('AA 滤波 + 抽取频谱 (f_s'' = %.1f MHz) - 无混叠', fs_dec/1e6));
grid on; xlim([-6 6]);
yl = ylim;

xline(-fs_dec/2/1e6, 'm--', 'LineWidth', 2);
xline( fs_dec/2/1e6, 'm--', 'LineWidth', 2);
text(fs_dec/2/1e6-0.5, yl(2)-5, sprintf('奈奎斯特\n边界'), ...
    'Color', 'm', 'FontSize', 9, 'FontWeight', 'bold');

% 标注保留的频率分量
text(0, yl(2)-15, sprintf('仅保留\n%.0f, %.0f, %.0f MHz', f1/1e6, f2/1e6, f3/1e6), ...
    'Color', 'g', 'FontWeight', 'bold', 'FontSize', 9, ...
    'HorizontalAlignment', 'center', 'BackgroundColor', [0.9 1 0.9]);
hold off;

%% ==================== 图 3: 频谱对比（关键对比图）====================
figure('Name', 'Key Comparison: Aliasing Effect');

subplot(1,2,1);
plot(f_dec, spec_dec, 'r', 'LineWidth', 2); hold on;
xline(-fs_dec/2/1e6, 'k--', 'LineWidth', 1.5);
xline( fs_dec/2/1e6, 'k--', 'LineWidth', 1.5);
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title('直接抽取 - 出现混叠');
grid on; xlim([-6 6]); ylim([-80 0]);
% 标注原始频率和混叠频率
text(3, -10, sprintf('3 MHz\n(混叠)'), 'Color', 'r', 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'center', 'BackgroundColor', [1 0.9 0.9]);
text(4, -20, '4 MHz', 'Color', 'b', 'FontWeight', 'bold');
legend('抽取后频谱', '奈奎斯特边界', 'Location', 'southwest');
hold off;

subplot(1,2,2);
plot(f_dec, spec_aa_dec, 'g', 'LineWidth', 2); hold on;
xline(-fs_dec/2/1e6, 'k--', 'LineWidth', 1.5);
xline( fs_dec/2/1e6, 'k--', 'LineWidth', 1.5);
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title('AA 滤波 + 抽取 - 无混叠');
grid on; xlim([-6 6]); ylim([-80 0]);
text(2, -10, '2 MHz', 'Color', 'g', 'FontWeight', 'bold');
text(4, -20, '4 MHz', 'Color', 'g', 'FontWeight', 'bold');
legend('抽取后频谱', '奈奎斯特边界', 'Location', 'southwest');
hold off;

%% ==================== 进阶：多相滤波实现（可选）====================
fprintf('\n进阶：多相滤波抽取验证\n');
baseband_tx_poly = upfirdn(baseband_tx, h_aa, 1, dsmp_rate);
% 补偿延迟并截取相同长度
baseband_tx_poly = baseband_tx_poly(delay_aa/dsmp_rate+1 : delay_aa/dsmp_rate+length(baseband_tx_aa_dec));
err = max(abs(baseband_tx_poly - baseband_tx_aa_dec));
fprintf('多相滤波与标准方法的最大误差: %.2e (应接近 0)\n', err);

多相滤波器

多相滤波器常用于这种多速率数字信号处理系统。

假设要对输入数据序列进行倍下采样，为了避免频谱展宽导致的混叠，首先要进行抗混叠滤波，然后对滤波结果进行下采样，显然，滤波输出中每个数据就有个数据会被丢弃，也就是说，针对这个数据的滤波处理是浪费的，多相滤波器就是为了解决这一计算浪费问题而设计的。

假设要对输入数据进行倍上采样，在上采样之后，需要进行滤波处理以消除镜像，一般做法是先进行插零上采样，然后进行滤波。采用多相滤波器可以颠倒采样和滤波的顺序以降低运算量。