Inverse filter (decorrelation)

Summary

For a minimum phase system, the exact inverse filter is generally an IIR filter.
It can be approximated with an FIR filter. This approach may also be used on non-minimum phase systems, if the compensation filter is allowed to delay the signal.

A matlab / Octave example can be found below.

Introduction

Frequently, the effect of one filter

should be compensated using another filter (Figure 1).

Figure 1: A filter compensating the frequency response of another filter

For example, a radio frequency transceiver or a measurement instrument may include analog filters that introduce significant passband distortion.

Exact inverse

A possible approach is to calculate the compensation filter as the exact inverse of

In many cases, the inverse filter is not a practical engineering solution.
It is still an interesting study case.

In the discrete-time domain, $H_{d}(f)$ becomes $H_{d}(z)$ , and

turns into

.
The transfer function of

can be written as:

$H_d(z)=\frac{b_0+b_{1}z^{-1}+b_{2}z^{-2}+\ldots}{1-a_{1}z^{-1}-a_{2}z^{-2}-\ldots}$

To compensate the frequency response of $H_{d}(z)$ , the product of

and

should be unity:

It follows that

$C(z)=\frac{1}{H_{d}(z)}$

In other words, the numerator of $H_{d}(z)$ becomes the denominator of

and vice versa.

Minimum phase requirement

A filter is minimum phase, if all its zeros fall inside the unit circle.
Inverting it turns zeros into poles.
Zeroes outside the unit circle become instable poles.
Therefore, the exact inverse of a non-minimum phase filter is always instable.

Nonzero amplitude response

Also, there may not be any zeros on the unit circle.
It is straightforward to see that once $H_{d}(z)$ becomes zero at a particular frequency, the signal is gone and cannot be recovered.
Note that such is the case for both ideal lowpass and highpass filters!
In other words: The amplitude response may not be zero anywhere “on the frequency axis”.

Approximate FIR inverse

Filters in technical applications are rarely ideal, and it may be sufficient to find an FIR filter that approximates the inverse to a given accuracy.
The matlab example program given below calculates the optimum compensation filter. It requires

the impulse response of $H_{d}(z)$
the nominal delay
the number of taps for the compensation filter

It will then find an equalizer impulse response as a least-squares optimum.

Inverting a filter with zeros on the frequency axis is dividing by zero.
Garbage in, garbage out.

Nominal delay

For a minimum phase system, the optimum delay is “none”:
Figure 2 shows the equalizer impulse response for such a case.
Nonetheless, the nominal delay was set to 6 samples. The least-squares solution gives zero for the first five taps before the main pulse, rendering them useless for equalization.

Figure 2: Impulse response of compensation filter for minimum-phase system

For a non-minimum phase system, the optimum delay should be “as long as possible”, shown in figure 3.
Here, the nominal delay is 45 samples. In the same manner as before, samples after the main pulse are not useful for the compensation filter:

Figure 3: Impulse response of compensation filter for non minimum-phase system

Finally, a system with zeros inside and outside the unit circle results in a compensation filter with an impulse response as shown in figure 4. Its envelope grows before the main tap, and decays afterwards.
The optimum delay depends on the time constants of the envelope.

Figure 4: Impulse response of compensation filter for zeros inside and outside the unit circle

Figure 5 shows the impulse response of the channel. The result after equalization is shown in figure 6, and again in figure 7 on a logarithmic scale as $20\log_{10}(c(t))$ .

Figure 5: Unequalized impulse response of the “channel”

Figure 6: Combined impulse response “channel” and equalizer

Figure 7: Combined impulse response “channel” and equalizer (logarithmic scale)

The arrow in figure 7 points to the wanted impulse (unity = 0 dB).

Code (Matlab / Octave)

Download

% *************************************************************
% Frequency response compensation
% © Markus Nentwig 2007
% provided without any warranty
%
% Calculation of a compensation filter to equalize the
% frequency response
% Note that this approach is not too useful for most
% engineering applications
% *************************************************************
close all; clear all;

% *************************************************************
% nomDelay: Latency
% --------
% if the channel is not minimum phase, the compensation filter
% becomes noncausal. This is taken into account by adding latency
% to the system.
% For a minimum phase channel, the compensation filter could be
% described by a finite number of coefficients. 
% The impulse response is causal, but still infinitely long.
%
% Number of taps
% --------------
% The impulse response is infinitely long in every non-trivial
% case. An implementation as FIR filter needs to truncate it.
% Therefore, the resulting compensation filter is only an
% approximation anyway.
% *************************************************************

if 0
  % Use conjugate pole pairs for no particular reason
  % zeros inside the unit circle - minimum phase system
  chan = poly([0.8*exp(0.3i*pi), 0.8*exp(-0.3i*pi)]);
  nomDelay=0; % optimum - all taps left of the main pulse remain zero
  
elseif 0
  % zeros outside the unit circle - nonminimum phase system
  chan = poly([1.2*exp(0.3i*pi), 1.2*exp(-0.3i*pi)]);
  nomDelay=1; % optimum - all taps right of the main pulse remain zero
  
elseif 1
  % zeros both inside and outside the unit circle- nonminimum phase system
  chan = poly([1.2*exp(0.3i*pi), 1.2*exp(-0.3i*pi), 0.8*exp(0.2i*pi), 0.8*exp(-0.2i*pi)]);
  % compensation filter impulse response is both left and right of the main pulse
  nomDelay=0.5; % arbitrary choice, not optimum
end

ntaps=50; % compensation filter length. Size -does- matter :)
% The convolution of channel and compensation filter should result in a pulse
% nomDelay gives the location of the pulse :
% 0: no delay (causal)
% 1: maximum delay (non-causal)
pulsePos=floor(nomDelay*(ntaps+length(chan)-1))+1;

% *************************************************************
% Calculate equalizer
% *************************************************************
base=[];
pchan=[chan zeros(1, ntaps)];
for k=1:ntaps
  base=[base; circshift(pchan, [0, k-1])];
end

% the desired impulse response of channel -and- equalizer
ideal=zeros(1, length(chan)+ntaps); 
ideal(pulsePos)=1; 

% Backslash operator finds the least-squares solution for:
% sum of
% - 1st EQ tap times undelayed channel impulse response
% - plus 2nd EQ tap times 1-delayed ch. IR
% - plus 3rd EQ tap times 2-delayed ch. IR
% - plus ...
% equals ideal impulse response, allowing for some delay
eq=transpose(base)\transpose(ideal);

% *************************************************************
% Test it and plot
% *************************************************************
v=[1]; % test pulse
%v=[1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6i]; % alternative test pulse

% ... through channel (convolve with impulse response)
v=conv(v, chan);

% ... through equalizer (convolve with impulse response)
v=conv(v, eq);

figure(); stem(abs(v));
title('equalized pulse (linear scale)');
figure(); plot(20*log10(abs(v)+1e-20), '+');
title('equalized pulse (logarithmic scale)'); ylabel('dB');
figure(); stem(eq);
title('equalizer impulse response');

Comments? Questions?

Please send me a mail! mnentwig@elisanet.fi