FIR滤波去除群时延的方法

志成于冰

已于 2024-06-19 09:26:23 修改

阅读量1.2k

点赞数 4

文章标签：信号处理 matlab 笔记

于 2024-06-18 17:39:40 首次发布

本文链接：https://blog.csdn.net/weixin_43853598/article/details/139778672

版权

什么是群时延

群时延是指信号在通过滤波器或传输系统时，不同频率成分的相位延迟变化率。
IIR（无限脉冲响应）滤波器使用反馈结构，其当前输出依赖于当前和过去的输入值以及过去的输出值。通常具有非线性的相位响应，不同频率分量会经历不同的相位延迟，导致群时延不恒定。
FIR（有限脉冲响应）滤波器使用非反馈结构，其输出仅依赖于当前和过去的输入值。FIR滤波器可以设计成具有线性相位响应，所有频率分量经历相同的相位延迟。如下图，滤波处理后出现时域偏移。

为什么要去除群时延

信号完整性：如果不同频率成分的延迟不同，信号的时域波形会被拉伸或压缩，导致信号失真。这对于需要高保真度的信号处理非常不利，例如音频处理、图像处理等。
相位失真：群时延变化会引起相位失真。例如，在相位调制的通信系统中，相位失真会直接影响解调结果，导致通信质量下降。
同步问题：在多通道信号处理中，不同通道的群时延差异会导致信号不同步，影响整体系统性能。例如在多麦克风阵列中，各麦克风的信号需要准确同步，否则会影响波束形成和方向定位的效果。

如何去除群时延

直接截取对齐：对信号进行切割，保证滤波前后的信号对齐。方法简单，但可能会导致信号丢失部分信息。
conv 卷积法：调用conv函数进行信号和滤波器参数的卷积运算。简单易用、计算速度快，但会引入滤波器本身的群时延，会在信号的边界处产生畸变。
filtfilt 双向滤波法：filtfilt函数通过前向和后向滤波，可以完全消除滤波器的群时延，相位保真度高。但是计算复杂度较高，是单向滤波计算量的两倍；对于较短的信号，可能会出现边界效应。
原始信号补零：先在原始信号后补零，再进行滤波处理。简单灵活、减小滤波器引入的群时延。但是补零后的信号在滤波过程中可能会引入一定的畸变，特别是在信号的边界处；计算复杂度增加。
滤波性能参数：
信号幅度保真度：滤波后信号的幅度变化。
相位保真度：滤波后信号的相位变化。
信号能量损失：滤波后信号的总能量变化。

%% Clear
clc; clear; close all;
%% 生成测试信号
fs = 125; % 采样频率
t = 0:1/fs:85-1/fs; % 时间向量
f1 = 1; % 频率成分1
f2 = 5; % 频率成分2
f3 = 10; % 频率成分2
signal = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t); % 生成信号
% 设计一个FIR滤波器
fcuts = [5  10];
mags = [1 0];
devs = [0.00001 0.05];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
b = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); 

%% 计算群时延
grpdelay(b,size(signal, 2), fs)
if(mod(n, 2) == 0)
    group_delay = n / 2;
else
    group_delay = (n-1) / 2;
end

%% filter滤波信号
filter_signal = filter(b, 1, signal);
figure;
plot(signal, '-r');
hold on
plot(filter_signal, '-b');
title('原始信号/filter滤波');
legend('raw','filter滤波');
xlabel('时间 (秒)');
ylabel('幅值');


%% 1. conv滤波信号，卷积会在信号的边界处产生效应
conv_signal = conv(signal, b, 'same');
figure;
plot(signal, '-r');
hold on
plot(conv_signal, '-b');
title('原始信号/conv滤波');
legend('raw','conv滤波');
xlabel('时间 (秒)');
ylabel('幅值');

%% 2. filtfilt滤波信号，无失真的滤波结果，双向滤波的计算量大约是其两倍
filtfilt_signal = filtfilt(b, 1, signal);
figure;
plot(signal, '-r');
hold on
plot(filtfilt_signal, '-b');
title('原始信号/filtfilt滤波');
legend('raw','filtfilt滤波');
xlabel('时间 (秒)');
ylabel('幅值');


%% 3. 补偿群时延--在信号后填充零点，计算复杂度增加,边界处可能会引入一定的畸变
new_signal = [signal, zeros(1, group_delay)];
padding_zeros_signal = filter(b, 1, new_signal);
padding_zeros_signal = padding_zeros_signal(group_delay+1:end);
figure;
plot(signal, '-r');
hold on
plot(padding_zeros_signal, '-b');
title('原始信号/filter滤波(原信号补零)');
legend('raw','filter滤波(原信号补零)');
xlabel('时间 (秒)');
ylabel('幅值');

%% 4. 补偿群时延--截取信号对齐，直接简单，丢失信号信息，如数据量较少可直接pass
figure;
plot(signal(1:end-group_delay), '-r');
hold on
plot(filter_signal(group_delay+1:end), '-b');
title('原始信号/filter滤波(截取信号对齐)');
legend('raw','filter滤波截取信号');
xlabel('时间 (秒)');
ylabel('幅值');

signal_para.raw_signal = signal;
signal_para.conv_signal = conv_signal;
signal_para.filtfilt_signal = filtfilt_signal;
signal_para.padding_zeros_signal = padding_zeros_signal;
compare_filter_performance(signal_para);

function [] = compare_filter_performance(signal_para)
    % 计算幅度最大误差
    amplitude_error_conv = max(abs(signal_para.raw_signal - signal_para.conv_signal));
    amplitude_error_filtfilt = max(abs(signal_para.raw_signal - signal_para.filtfilt_signal));
    amplitude_error_zero_padded = max(abs(signal_para.raw_signal - signal_para.padding_zeros_signal));
    amplitude_error = [amplitude_error_conv, amplitude_error_filtfilt, amplitude_error_zero_padded];

    % 计算能量损失
    energy_loss_conv = (sum(signal_para.raw_signal.^2) - sum( signal_para.conv_signal.^2))/size(signal_para.raw_signal, 2);
    energy_loss_filtfilt = (sum(signal_para.raw_signal.^2) - sum(signal_para.filtfilt_signal.^2))/size(signal_para.raw_signal, 2);
    energy_loss_zero_padded = (sum(signal_para.raw_signal.^2) - sum(signal_para.padding_zeros_signal.^2))/size(signal_para.raw_signal, 2);    
    energy_loss = [energy_loss_conv,energy_loss_filtfilt,energy_loss_zero_padded];

    % 比较相位平均误差
    % 计算相位谱
    raw_phase = angle(fft(signal_para.raw_signal));
    conv_filtered_phase = angle(fft(signal_para.conv_signal));
    filtfilt_filtered_phase = angle(fft(signal_para.filtfilt_signal));
    zero_padded_filtered_phase = angle(fft(signal_para.padding_zeros_signal));
    % 计算相位平均误差
    phase_diff = [sum(abs(conv_filtered_phase-raw_phase)), sum(abs(filtfilt_filtered_phase - raw_phase)), sum(abs(zero_padded_filtered_phase-raw_phase))];
    phase_diff = phase_diff/size(raw_phase, 2);

    % 绘图对比幅度最大误差和能量损失
    data = [amplitude_error;energy_loss;phase_diff];
    plot_bar(data)    
end


function [] = plot_bar(signal)
    figure;
    hold on;
    box on;
    xLables = ["幅度最大误差", "能量平均损失","相位平均误差"];
    xGroup = ["conv处理", "filtfilt处理", "filter处理(原信号补零)"];
         
    bar(signal ,'FaceColor','flat');
    set(gca,'XTick', [1,2,3] , 'XTickLabel',xLables);
    xlabel("性能参数",'FontName','宋体','FontSize',14,'Fontweight','bold');                        
    hl = legend(xGroup);
    set(hl,'Box','off','location','NorthWest','NumColumns',1,'FontSize',12,'FontName','宋体');    
    ylabel('差值');
    title('三种滤波方法性能对比');
end