数字信号处理系列:多速率内插过程中镜像频率如何处理?
本节目录
一、多速率内插原理
1、时域分析
2、频域分析
3、低通滤波器抑制镜像
二、Matlab设计验证
三、Matlab仿真结果
本节内容
一、多速率内插原理
多速率内插(升采样)是信号处理中提高信号采样率的核心技术,通过插入额外样本点来增加信号的采样率,同时保持信号频谱的完整性。
1、时域分析
时域操作:零值内插(Zero-Padding)
原始信号为 x[n],采样率为fs,内插因子为 L(整数),则内插后的信号xup[n]:
时域波形
2、频域分析
频域复制效应,零值内插在频域表现为原始频谱的周期性重复。 原始信号x[n]的频谱为X(ejω),内插后的频谱Xup(ejω)为Xup(ejω)=X(ejωL)
混叠风险
3、低通滤波器抑制镜像
镜像抑制滤波器(Anti-Imagining Filter)
滤波过程
二、Matlab设计验证
参数设置
生成测试信号
内插操作
设计镜像抑制滤波器
频谱分析函数
plot_spectrum函数用于绘制信号的频谱,并可以添加频率标记。
%% 频谱分析函数
function plot_spectrum(signal, fs, title_str, freq_markers, marker_labels)
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;
% 绘图
figure;
plot(f, P1);
title(title_str);
xlabel('频率 (Hz)');
ylabel('幅度');
grid on;
xlim([0 fs/2]);
% 添加频率标记
if nargin > 3
hold on;
colors = ['r', 'g', 'b', 'm', 'c'];
for i = 1:length(freq_markers)
xline(freq_markers(i), '--', marker_labels{i}, 'Color', colors(i), 'LineWidth', 1.5);
end
hold off;
end
end
matlab具体代码如下:
%% 信号内插与频谱特性仿真
clear; clc; close all;
%% 参数设置
fs = 1000; % 原始采样率 (Hz)
T = 1; % 信号时长 (秒)
L = 3; % 内插因子
% 信号参数
f1 = 50; % 低频分量
f2 = 150; % 高频分量
f3 = 300; % 更高频分量 (在原始Nyquist频率内)
%% 生成测试信号
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_upsampled = zeros(1, N*L);
signal_upsampled(1:L:end) = signal;
fs_upsampled = fs * L;
% 设计镜像抑制滤波器
f_cutoff = fs/2 * 0.9; % 截止频率 (90% 原始Nyquist频率)
filter_order = 100; % 滤波器阶数
% FIR滤波器设计
b = fir1(filter_order, f_cutoff/(fs_upsampled/2), 'low');
% 应用低通滤波器
filtered_upsampled = filtfilt(b, 1, signal_upsampled);
%% 关键频率定义
nyquist_original = fs/2;
nyquist_upsampled = fs_upsampled/2;
freq_markers = [f1, f2, f3, nyquist_original, nyquist_upsampled];
marker_labels = {'f1=50Hz', 'f2=150Hz', 'f3=300Hz', ...
sprintf('f_{Nyq,orig}=%dHz', nyquist_original), ...
sprintf('f_{Nyq,up}=%dHz', nyquist_upsampled)};
%% 频谱分析
% 原始信号频谱
plot_spectrum(signal, fs, sprintf('原始信号频谱 (f_s = %d Hz)', fs), ...
freq_markers, marker_labels);
% 内插零值后频谱
plot_spectrum(signal_upsampled, fs_upsampled, ...
sprintf('内插零值信号频谱 (L=%d, f_s=%d Hz)', L, fs_upsampled), ...
freq_markers, marker_labels);
% 滤波后内插信号频谱
plot_spectrum(filtered_upsampled, fs_upsampled, ...
sprintf('滤波后内插信号频谱 (L=%d, f_s=%d Hz)', L, fs_upsampled), ...
freq_markers, marker_labels);
%% 时域信号可视化
t_upsampled = 0:1/fs_upsampled:T-1/fs_upsampled;
figure('Position', [100, 100, 1200, 900]);
% 原始信号
subplot(4,1,1);
plot(t, signal);
title(sprintf('原始信号 (f_s = %d Hz)', fs));
xlabel('时间 (s)');
ylabel('幅度');
xlim([0 0.05]);
% 内插零值后信号
subplot(4,1,2);
stem(t_upsampled, signal_upsampled, 'Marker', 'none');
title(sprintf('内插零值信号 (L=%d, f_s=%d Hz)', L, fs_upsampled));
xlabel('时间 (s)');
ylabel('幅度');
xlim([0 0.05]);
% 滤波后内插信号
subplot(4,1,3);
plot(t_upsampled, filtered_upsampled);
title(sprintf('滤波后内插信号 (L=%d, f_s=%d Hz)', L, fs_upsampled));
xlabel('时间 (s)');
ylabel('幅度');
xlim([0 0.05]);
% 原始信号与滤波后内插信号对比
subplot(4,1,4);
plot(t, signal, 'b', 'LineWidth', 2); hold on;
% 由于内插后采样率提高,需要下采样以与原始信号对齐
plot(t, filtered_upsampled(1:L:length(filtered_upsampled)), 'r--', 'LineWidth', 1.5);
title('原始信号与滤波后内插信号对比');
xlabel('时间 (s)');
ylabel('幅度');
legend('原始信号', '滤波后内插信号下采样');
xlim([0 0.05]);
%% 滤波器响应分析
figure('Position', [100, 100, 1000, 400]);
freqz(b, 1, 1024, fs_upsampled);
title('镜像抑制滤波器频率响应');
grid on;
%% 镜像频率分析
% 关键频率
original_freqs = [f1, f2, f3];
fprintf('\n镜像频率分析:\n');
% 对于每个分量
for k = 1:length(original_freqs)
f_orig = original_freqs(k);
% 计算镜像位置 (k * fs ± f_orig)
images = [];
for n = 1:L-1
image1 = n * fs - f_orig;
image2 = n * fs + f_orig;
images = [images, image1, image2];
end
% 在频谱中确认
[P1, f] = calculate_spectrum(signal_upsampled, fs_upsampled);
fprintf('原始频率 %d Hz:\n', f_orig);
for img = images
if img > 0 && img < fs_upsampled/2
% 找到最接近的频率点
[~, idx] = min(abs(f - img));
magnitude = P1(idx);
fprintf(' 镜像频率 %d Hz 处幅度: %.4f\n', img, magnitude);
end
end
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
三、Matlab仿真结果
频谱分析
原始信号频谱
内插零值后信号频谱
滤波后内插信号频谱
时域信号分析
滤波器响应分析
信号内插过程中的频谱变化以及镜像抑制滤波器的作用,通过内插零值,采样率提高到原来的L倍,但引入了镜像频谱。通过设计合适的低通滤波器,可以抑制这些镜像分量,从而得到采样率提高且不失真的信号。在时域上,滤波过程将零值替换为插值样本,使信号波形变得平滑。