Inverse of filter function

How to write the inverse of the filter function?

I use a filter function with input x to get output y. I want to get input x back.

y = filter(b,a,x,zi)
y = <inverse filter to get original input x>

Answer 1

[Christoph F.]

In the z domain, the transfer function of a filter H(z) is B(z)/A(z).

The inverse of the transfer function is A(z)/B(z). However, this only makes sense if the zeros of B(z) are inside the unit circle; if they are outside, the inverse filter will have poles outside the unit circle and hence be unstable.

[Paul]

In addition to the possiblity of the inverse filter having poles outside the unit circle, which is of practical concern, another consideration is the causality of the inverse filter.

Define the signal

rng(100);
N = 1000;
sig = sin(2*pi*40*(0:N-1)*0.001);
n = 0:N-1;

Define an IIR filter with one zero and four poles and filter the signal

[B,A] = zp2tf(0.6,[0.1 0.2 0.3 0.4],1)

B = 1×5

0 0 0 1.0000 -0.6000

A = 1×5

1.0000 -1.0000 0.3500 -0.0500 0.0024

sigf  = filter(B,A,sig);

figure
plot(n,sig,n,sigf)

Attempting to filter with the inverse filter casues an error because the inverse filter is non-causal (more zeros than poles)

try
    sigr = filter(A,B,sigf);
catch ME
    ME.message
end

ans = 'First denominator filter coefficient must be non-zero.'

Instead, we mutiply the inverse filter by 1/z^3 (the power is the number of leading zeros in B), then filter

sigr = filter(A,B(end-1:end),sigf);

then shift the output back by three samples.

sigr = sigr(4:end);
nr   = numel(sigr);

This approach recovered the input signal except for the last three samples.

figure
plot(n,sig,n(1:nr),sigr)

figure
plot(n(1:nr),sig(1:nr)-sigr)

Answer 2

[Jan]

This is not possible. Remember that filtering destroys the original information. If you e.g. remove the high frequency noise from a signal, it does not exist in the output anymore. Then there is no way to re-create it.

[Christoph F.]

Except for zeros in the amplitude response, signal components do not get removed, merely attenuated.

Mathematically, it is possible to invert B(z)/A(z) transfer function as A(z)/B(z), with some precautions (like B(z) having zeros outside the unit circle resulting in an unstable inverted filter). Profound knowledge of the mathematical background of digital signal processing is extremely helpful here, and unfortunately such background cannot be compressed into a short posting. Information can be found e.g. in Proakis/Manolakis "Digital Signal Processing" in the sub-chapter "Invertibility of Linear Time-Invariant Systems".

[Jan]

But there is no procedure, which solves this in general:

y = filter(b,a,x,zi)
x = <inverse filter to get original input y>

b,a could define a moving average filter. Then a reconstruction of the original signal is not possible. Or do I miss a point here?

[Christoph F.]

Just like there is no general inverse operation to a multiplication, i.e.

y=b*x
x=F(y) for any given y including zero?

Anything that was multiplied by zero cannot be restored. That does not reduce the usefulness of division as the inverse operation to multiplication.

[John D'Errico]

Jan, consider this simple problem:

K = ones(1,3)/3;
S0 = rand(1,10);
S1 = conv(S0,K,'same');

Given only S1 and the kernel K, can you recover the original signal? (To within floating point precision.) I've written it in the form of a conv question, but since filter and conv are really the same thing under the hood, this still applies.

Lets try, and see what happens.

We can write the convolution in terms of a matrix multiply.

A = spdiags(repmat(1/3,10,3),-1:1,10,10);

A is full rank.

rank(full(A))
ans =
    10

S1m = S0*A;
S1 - S1m
ans =
   0     0     0     0     0     0     0     0     0     0

And as we see, we accomplish the same result by writing the convolution as a matrix multiplication. And as long as the convolution matrix is nonsingular, we can recover the input. Here, slash applies, since the problem has been posed in linear algebraic form.

S0m = S1/A;
S0 - S0m
ans =
  -2.2204e-16   2.2204e-16            0  -1.1102e-16   1.1102e-16  -2.2204e-16   1.1102e-16            0  -5.5511e-17   8.3267e-17

A is even pretty well conditioned.

cond(full(A))
ans =
     17.255

so we don't expect to lose much of the original signal. If it were true that A was not full rank, then the original signal is not recoverable. Or, if A has a large condition number, then there will be numerical issues. Both cases can be viewed from a signal processing point of view too.

[Christoph F.]

A short example:

% Generate test signal and filter
t=0:0.001:2;
sig=sin(2*pi*40*t);
[B, A]=ellip(3, 0.5, 60, 39/500);
% Plot test signal, filtered signal, recovered test signal
plot(t, sig, 'k-', t, filter(B, A, sig), 'b-', t, filter(A, B, filter(B, A, sig)), 'r--');
% Maximum difference between recovered and test signal:
max(abs(sig - filter(A, B, filter(B, A, sig))))

The maximum difference of test signal and recovered test signal is in the order of 10^-11 in this example. Be aware that A(z)/B(z) has poles on the unit circle and is therefore only marginally stable.

[Jan]

Thanky you, John and Christoph.

[John D'Errico]

I think the difference is, suppose you took the mean of your data. Then the original data is not recoverable. But here you have a linear transformation of the data. N points in, N points out. As long as the mapping is not singular, then there is an inverse transformation.

[Paul]

@John D'Errico

Using the 'same' option with conv is not the same operation as using filter. In this example, the 'same' option results in an output that is shifte by 1 sample relative to the output of filter. Can the A matrix be changed appropriately?

Probably won't change the general conclusion ....

rng(100)
K = ones(1,3)/3;
S0 = rand(1,10);
S1 = conv(S0,K,'same');
S1a = conv(S0,K);
S1b = filter(K,1,S0);
[S1;S1a(1:numel(S0));S1b]

ans = 3×10

0.2739 0.4154 0.5159 0.4247 0.3237 0.2657 0.5394 0.5444 0.5126 0.2373 0.1811 0.2739 0.4154 0.5159 0.4247 0.3237 0.2657 0.5394 0.5444 0.5126 0.1811 0.2739 0.4154 0.5159 0.4247 0.3237 0.2657 0.5394 0.5444 0.5126