数字信号处理系列:多速率抽取过程中抗混叠滤波器设计
本节目录
一、多速率抽取原理
二、Matlab设计验证
三、Matlab仿真结果
本节内容
一、多速率抽取原理
多速率信号处理中的抽取（Decimation）是一种降低信号采样率的技术，其核心原理是通过有选择地保留部分样本来减少数据量，同时避免信息损失和频谱混叠。
抽取的基本原理
抽取操作可以表示为y[n]=x[D·n]，其中，D是抽取因子（整数），表示每隔D个样本取一个样本。
为什么需要抗混叠滤波？
直接抽取可能导致频谱混叠（Aliasing）。根据奈奎斯特采样定理，如果原始信号中包含高于新采样率一半的频率分量（即高于fs/(2D)），则这些高频分量在抽取后会混叠到低频区域，造成失真。
为了避免混叠，需要在抽取前进行低通滤波（抗混叠滤波），滤除高于新Nyquist频率（fs/(2D)）的分量。原始信号频谱周期为fs，抽取后频谱周期变为fs/D，若不滤波，高于fs/(2D)的分量将混叠到[0,fs/(2D)]区间。
抗混叠滤波器设计
二、Matlab设计验证
通过MATLAB仿真，验证信号在抽取过程中频谱混叠与不混叠的现象，并分析抗混叠滤波器的效果。
参数设置：
信号生成：
信号抽取与抗混叠滤波
matlab具体代码如下：
%% 信号抽取与频谱混叠仿真验证
clear; clc; close all;
%% 参数设置
fs = 1000; % 原始采样率 (Hz)
T = 1; % 信号时长 (秒)
D = 4; % 抽取因子
% 信号参数
f1 = 50; % 低频分量 (Hz) - 不会混叠
f2 = 150; % 高频分量 (Hz) - 会混叠
f3 = 300; % 更高频分量 (Hz) - 严重混叠
%% 生成测试信号
t = 0:1/fs:T-1/fs; % 时间向量
N = length(t); % 样本点数
% 多分量信号
signal = sin(2*pi*f1*t) + 0.5*sin(2*pi*f2*t) + 0.3*sin(2*pi*f3*t);
%% 抽取信号
% 直接抽取（无抗混叠滤波）
signal_decimated = signal(1:D:end);
fs_decimated = fs/D; % 新采样率
t_decimated = t(1:D:end);
%% 抗混叠滤波后抽取
% 设计抗混叠滤波器
f_cutoff = fs_decimated/2 * 0.9; % 截止频率 (90% Nyquist)
filter_order = 100; % 滤波器阶数
% FIR滤波器设计
b = fir1(filter_order, f_cutoff/(fs/2), 'low');
filtered_signal = filtfilt(b, 1, signal); % 零相位滤波
% 滤波后抽取
filtered_decimated = filtered_signal(1:D:end);
%% 频谱分析
% 原始信号频谱
plot_spectrum(signal, fs, sprintf('原始信号频谱 (f_s = %d Hz)', fs));
% 直接抽取信号频谱
plot_spectrum(signal_decimated, fs_decimated, ...
sprintf('直接抽取信号频谱 (D=%d, f_s=%d Hz)', D, fs_decimated));
% 抗混叠滤波后抽取信号频谱
plot_spectrum(filtered_decimated, fs_decimated, ...
sprintf('抗混叠滤波后抽取信号频谱 (D=%d, f_s=%d Hz)', D, fs_decimated));
%% 混叠现象可视化
% 计算Nyquist频率
nyquist_original = fs/2;
nyquist_decimated = fs_decimated/2;
% 创建频率标记
freq_markers = [f1, f2, f3, nyquist_decimated, nyquist_original];
marker_labels = {'f1=50Hz', 'f2=150Hz', 'f3=300Hz', ...
sprintf('f_{Nyq,dec}=%dHz', nyquist_decimated), ...
sprintf('f_{Nyq,orig}=%dHz', nyquist_original)};
% 绘制所有频谱对比
figure('Position', [100, 100, 1200, 800]);
% 原始信号频谱
subplot(3,1,1);
[P1, f] = calculate_spectrum(signal, fs);
plot(f, P1);
title(sprintf('原始信号频谱 (f_s = %d Hz)', fs));
xlabel('频率 (Hz)');
ylabel('幅度');
grid on;
xlim([0 fs/2]);
add_freq_markers(freq_markers, marker_labels);
% 直接抽取信号频谱
subplot(3,1,2);
[P1_dec, f_dec] = calculate_spectrum(signal_decimated, fs_decimated);
plot(f_dec, P1_dec);
title(sprintf('直接抽取信号频谱 (D=%d, f_s=%d Hz)', D, fs_decimated));
xlabel('频率 (Hz)');
ylabel('幅度');
grid on;
xlim([0 fs_decimated/2]);
add_freq_markers(freq_markers, marker_labels);
% 抗混叠滤波后抽取信号频谱
subplot(3,1,3);
[P1_filt, f_filt] = calculate_spectrum(filtered_decimated, fs_decimated);
plot(f_filt, P1_filt);
title(sprintf('抗混叠滤波后抽取信号频谱 (D=%d, f_s=%d Hz)', D, fs_decimated));
xlabel('频率 (Hz)');
ylabel('幅度');
grid on;
xlim([0 fs_decimated/2]);
add_freq_markers(freq_markers, marker_labels);
%% 混叠频率计算验证
% 对于f2=150Hz分量
f2_aliased = abs(f2 - fs_decimated * round(f2/fs_decimated));
fprintf('\n高频分量混叠验证:\n');
fprintf('原始频率 f2 = %d Hz\n', f2);
fprintf('抽取后理论混叠频率 = %d Hz\n', f2_aliased);
% 在频谱中找到实际混叠频率
[~, idx] = max(P1_dec);
f_measured = f_dec(idx);
fprintf('频谱中最高峰频率 = %.2f Hz\n', f_measured);
%% 时域信号可视化
figure('Position', [100, 100, 1200, 600]);
% 原始信号
subplot(3,1,1);
plot(t, signal);
title('原始信号');
xlabel('时间 (s)');
ylabel('幅度');
xlim([0 0.1]); % 显示前0.1秒
% 直接抽取信号
subplot(3,1,2);
plot(t_decimated, signal_decimated, 'o-');
title('直接抽取信号');
xlabel('时间 (s)');
ylabel('幅度');
xlim([0 0.1]);
% 抗混叠滤波后抽取信号
subplot(3,1,3);
plot(t_decimated, filtered_decimated, 'o-');
title('抗混叠滤波后抽取信号');
xlabel('时间 (s)');
ylabel('幅度');
xlim([0 0.1]);
%% 滤波器响应分析
figure('Position', [100, 100, 1000, 400]);
freqz(b, 1, 1024, fs);
title('抗混叠滤波器频率响应');
grid on;
频谱绘图函数
%% 频谱分析函数
function plot_spectrum(signal, fs, title_str)
N = length(signal);
NFFT = 2^nextpow2(N); % 最接近的2的幂次
% 计算频谱
Y = fft(signal, NFFT);
P2 = abs(Y/NFFT);
P1 = P2(1:NFFT/2+1);
P1(2:end-1) = 2*P1(2:end-1); % 单边谱
% 频率向量
f = fs*(0:(NFFT/2))/NFFT;
% 绘图
figure;
plot(f, P1);
title(title_str);
xlabel('频率 (Hz)');
ylabel('幅度');
grid on;
xlim([0 fs/2]);
end
频谱分析函数
function [P1, f] = calculate_spectrum(signal, fs)
N = length(signal);
NFFT = 2^nextpow2(N);
Y = fft(signal, NFFT);
P2 = abs(Y/NFFT);
P1 = P2(1:NFFT/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(NFFT/2))/NFFT;
end
辅助标记函数
function add_freq_markers(freqs, labels)
hold on;
colors = ['r', 'g', 'b', 'm', 'c'];
for i = 1:length(freqs)
xline(freqs(i), '--', labels{i}, 'Color', colors(i), 'LineWidth', 1.5);
end
hold off;
end
三、Matlab仿真结果
通过plot_spectrum函数对信号进行FFT分析，绘制频谱图。