数字信号处理中的逆滤波器（inverse filter）matlab实现

由线性时不变系统的可逆性，输入经过原系统之后再经过对应的逆系统可得到原输入。

用冲击响应来描述原系统和逆系统的关系，设原系统冲击响应为, 逆系统的冲激响应为，上图等价于如下的卷积方程：

这意味着：

上式可求出给定的，然而在时域中求解上式通常很困难。更简单的方法是将上式转换到域来求解逆系统，上式在变换域可表示为：

逆系统函数为：

的所有零点都位于单位圆内（最小相位系统），这样得到的是稳定的。

知道求似乎很简单，但是如果是一个很长的脉冲响应（FIR滤波器的系数），对其求逆系统会得到全极点的IIR滤波器，并且阶数也非常大。

在MIT的HRTF 数据库中，为了进行扬声器均衡，作者用到了使用DFT/IDFT进行求解扬声器脉冲响应的逆滤波器的方法。具体的步骤如下：

1. 对进行padding（补零），使用FFT计算的离散傅里叶变换（DFT）得到。补零是防止发生混叠；
2. 计算，限制幅度
3. 使用IFFT计算的IDFT
4. 得到因果的滤波器，并加窗

matlab代码如下

[invdft, fs3] = audioread('.\data\headphones+spkr\Opti-inverse.wav');
[spkir, fs4] = audioread('.\data\headphones+spkr\Optimus.wav');


N = 16384;
Y = fft(spkir, N);
ind1 = floor(18000/44100*N);
mag = abs(Y);
ang = angle(Y);

% mag(1) = 1 / mag(1);
mag(2:ind1) = 1 ./ mag(2:ind1);
mag(ind1+1:N-ind1+1) = mag(ind1);
mag(N-ind1+2:N) = flipud(mag(2:ind1));

% mag(mag2db(mag) > 60) = db2mag(60);

Ynew1 = mag.*exp(1j*ang);

yh = fftshift(real(ifft(Ynew1)));
[~, loc] = max(yh);
yhs = yh(loc-1023:loc+1024).*hann(2048)*db2mag(-6);
figure;
plot(yhs);
title('inverse filter')
figure;
freqz(yhs,1);
title('inverse filter frequency response')

fvtool(yhs, 1, invdft,1)

Optimus.wav是MIT HRTF数据集中的扬声器脉冲响应，长度为16383。Opti-inverse.wav是MIT数据库提供的2048点的扬声器脉冲响应的逆滤波器系数，上面的代码以Optimus.wav作为输入，求解逆滤波器。最终得到的逆滤波器和数据库提供的比较接近 ，如下图：

参考：