Out of curiosity, is it possible to actually reverse a filtering on MATLAB. For example, if I have

[b,a] = butter(6,0.6);
dataIn = randn(1000,1);
dataOut = filter(b,a,dataIn);

Is there some way to go from dataOut and obtain dataIn in this example? I was able to more or less get there with the very simple

dataOutprime = filter(a,b,dataOut);

dataOutprime still doesn't exactly overlay on top of dataIn though.


In general it's only possible to implement causal and stable filters. There are exceptions where marginally stable filters are used, but this doesn't apply here. So if you want to invert a given filter, this is only possible if its zeros are inside the unit circle of the complex plane (such filters are called minimum-phase filters). The zeros of the given filter become the poles of the inverse filter, so if the zeros are inside the unit circle, the poles of the inverse filter will also be inside the unit circle, and, consequently, the resulting filter is causal and stable.

The problem with the filter in your example is that it has a sixth order zero at $z=-1$, i.e. on the unit circle. Consequently, this filter is not strictly minimum-phase and it can't be inverted by a causal and stable filter. However, it can be inverted approximately by a causal and stable filter. One way to achieve such an approximate inversion is by minimizing the mean squared error between the ideal and the actual filter response. In practice, it is often acceptable to add some delay to the output, which will make it easier to approximate the desired response by a causal and stable filter.

Note that transients don't play a role. It's the inverse filter's instability that is the problem here.

As an example, take a filter with coefficients

b = [1,0,-0.5]; a = [1,0.3,0.1];

The zeros of the numerator polynomial (with coefficients b) are all inside the unit circle, so the inverse filter is causal and stable.

x = rand(100,1);
y = filter(b,a,x);
x2 = filter(a,b,y);
max(abs(x2-x))    % in the order of 1e-16, i.e., machine precision
| improve this answer | |

I think the general approach is correct. You don't get the same result due to numerical instability of high order IIR filters. You may try breaking filter into cascade of 2 order filters and do the same afterwards (swapping a & b)

| improve this answer | |
  • $\begingroup$ It is also due to the transient response of the filters. You can improve on the numerical "instability" (quantization error) by switching to a larger wordsize, and using a balanced realization or lattice ladder implementation. $\endgroup$ – Arnfinn Jun 24 '16 at 7:48
  • 1
    $\begingroup$ @Arnfinn: Transients are no problem; if you can invert the filter by a causal and stable filter, then the inversion is perfect. $\endgroup$ – Matt L. Jun 24 '16 at 10:03
  • $\begingroup$ @VK69UK: The problem here is the fact that there exists no causal and stable filter that exactly inverts the given filter. It has nothing to do with numerical instabilities of high order IIR filters, so breaking it up in 2nd order sections won't help. $\endgroup$ – Matt L. Jun 24 '16 at 10:09
  • $\begingroup$ @MattL.: Sampling zeros, that's right... Use this: researchgate.net/publication/… $\endgroup$ – Arnfinn Jun 24 '16 at 11:55
  • $\begingroup$ @Arnfinn: Sampling zeros? $H(z)/H(z)=1$ (if $1/H(z)$ can be implemented). $\endgroup$ – Matt L. Jun 24 '16 at 13:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.