Inverse of a simple non-linear-phase FIR filter

Hello everyone,

   Can anyone tell me if there is a known way to find the inverse of a
given FIR filter? The idea is that given some signal y(x), and an FIR
filter f(x), I want to
find a way to compute g(x) which is the inverse of f(x) such that:

      y(x) * f(x) * g(x) = y(x).  
      (* = convolution)

Meaning if I have some signal, I can apply an FIR filter to shape it
then apply an inverse FIR to bring back the original signal.

Any ideas?
--Shafik

shafik@u.arizona.edu (Shafik) wrote in message news:<5e5cc96b.0408182053.6c5bab12@posting.google.com>...
> Hello everyone,
> 
>    Can anyone tell me if there is a known way to find the inverse of a
> given FIR filter? The idea is that given some signal y(x), and an FIR
> filter f(x), I want to
> find a way to compute g(x) which is the inverse of f(x) such that:
> 
>       y(x) * f(x) * g(x) = y(x).  
>       (* = convolution)
> 
> Meaning if I have some signal, I can apply an FIR filter to shape it
> then apply an inverse FIR to bring back the original signal.
> 
> Any ideas?

You can get what you want, provided the FIR filter is minimum phase, i.e.
has all its zeros in the interior of the unit circle. In that case, you 
can get a stable, causal IIR filter that compensates the effects of 
your FIR filter. 

See e.g. section 4.6 in

Proakis & Manolakis: "Digital Signal processing. Pronciples, Algorithms
   and Applications" Prentice-Hall, 3rd ed. 1996.

Rune

Rune,

As long as there are no zeros on the 
unit circle, it seems to me the FIR
filter is invertible, even if it 
is not minimum phase.  Minimum phase
means a causal inverse, but if the 
entire signal can be buffered before processing,
then it seems to me that the inverse does not
need to be causal. Am I missing something?

In article <f56893ae.0408190125.27ec4942@posting.google.com>,
Rune Allnor <allnor@tele.ntnu.no> wrote:
>shafik@u.arizona.edu (Shafik) wrote in message news:<5e5cc96b.0408182053.6c5bab12@posting.google.com>...
>> Hello everyone,
>> 
>>    Can anyone tell me if there is a known way to find the inverse of a
>> given FIR filter? The idea is that given some signal y(x), and an FIR
>> filter f(x), I want to
>> find a way to compute g(x) which is the inverse of f(x) such that:
>> 
>>       y(x) * f(x) * g(x) = y(x).  
>>       (* = convolution)
>> 
>> Meaning if I have some signal, I can apply an FIR filter to shape it
>> then apply an inverse FIR to bring back the original signal.
>> 
>> Any ideas?
>
>You can get what you want, provided the FIR filter is minimum phase, i.e.
>has all its zeros in the interior of the unit circle. In that case, you 
>can get a stable, causal IIR filter that compensates the effects of 
>your FIR filter. 
>
>See e.g. section 4.6 in
>
>Proakis & Manolakis: "Digital Signal processing. Pronciples, Algorithms
>   and Applications" Prentice-Hall, 3rd ed. 1996.
>
>Rune

Is there any way you can show me a simple example? Im new to this :-(

Thanks
--Shafik



gelbart@ICSI.Berkeley.EDU (David Gelbart) wrote in message news:<cg3069$gvm$1@agate.berkeley.edu>...
> Rune,
> 
> As long as there are no zeros on the 
> unit circle, it seems to me the FIR
> filter is invertible, even if it 
> is not minimum phase.  Minimum phase
> means a causal inverse, but if the 
> entire signal can be buffered before processing,
> then it seems to me that the inverse does not
> need to be causal. Am I missing something?
> 
> In article <f56893ae.0408190125.27ec4942@posting.google.com>,
> Rune Allnor <allnor@tele.ntnu.no> wrote:
> >shafik@u.arizona.edu (Shafik) wrote in message news:<5e5cc96b.0408182053.6c5bab12@posting.google.com>...
> >> Hello everyone,
> >> 
> >>    Can anyone tell me if there is a known way to find the inverse of a
> >> given FIR filter? The idea is that given some signal y(x), and an FIR
> >> filter f(x), I want to
> >> find a way to compute g(x) which is the inverse of f(x) such that:
> >> 
> >>       y(x) * f(x) * g(x) = y(x).  
> >>       (* = convolution)
> >> 
> >> Meaning if I have some signal, I can apply an FIR filter to shape it
> >> then apply an inverse FIR to bring back the original signal.
> >> 
> >> Any ideas?
> >
> >You can get what you want, provided the FIR filter is minimum phase, i.e.
> >has all its zeros in the interior of the unit circle. In that case, you 
> >can get a stable, causal IIR filter that compensates the effects of 
> >your FIR filter. 
> >
> >See e.g. section 4.6 in
> >
> >Proakis & Manolakis: "Digital Signal processing. Pronciples, Algorithms
> >   and Applications" Prentice-Hall, 3rd ed. 1996.
> >
> >Rune

"Shafik" <shafik@u.arizona.edu> wrote in message
news:5e5cc96b.0408191832.3060635b@posting.google.com...
> Is there any way you can show me a simple example? Im new to this :-(
>
> Thanks
> --Shafik

It's probably not very useful.  Unless the original FIR is minimum-phase,
the inverse filter will be an unstable IIR.  If you really want to do it,
then it's easy.  In DSP-talk, you just invert the transfer function.  A more
algorithmic description:

Set up the original FIR to process your *output* signal.  At every sample
instant, calculate the new output sample value that will produce a new
filtered sample that matches the corresponding input sample.  This is easily
done using the leading coefficient of the FIR, which must not be 0, and the
current FIR state.  Since you have ensured that filtering the output signal
with the original FIR would produce the input signal, you know you have
applied the FIR's inverse.

--
Matt

gelbart@ICSI.Berkeley.EDU (David Gelbart) wrote in message news:<cg3069$gvm$1@agate.berkeley.edu>...
> Rune,
> 
> As long as there are no zeros on the 
> unit circle, it seems to me the FIR
> filter is invertible, even if it 
> is not minimum phase.  Minimum phase
> means a causal inverse, but if the 
> entire signal can be buffered before processing,
> then it seems to me that the inverse does not
> need to be causal. Am I missing something?

You are right. I don't want to bring in buffering of off-line 
issues until they are specifically asked for, or I know enough 
about the application to know they are applicable. 

Rune

Hi Shafik,

My thinking relates to section 5.2 of Discrete-Time Signal Processing
(2nd edition) by Oppenheim and Schafer.  There are some examples in
chapter 5
that you may find helpful.   

shafik@u.arizona.edu (Shafik) wrote in message news:<5e5cc96b.0408191832.3060635b@posting.google.com>...
> Is there any way you can show me a simple example? Im new to this :-(
> 
> Thanks
> --Shafik
> 
> 
> 
> gelbart@ICSI.Berkeley.EDU (David Gelbart) wrote in message news:<cg3069$gvm$1@agate.berkeley.edu>...
> > Rune,
> > 
> > As long as there are no zeros on the 
> > unit circle, it seems to me the FIR
> > filter is invertible, even if it 
> > is not minimum phase.  Minimum phase
> > means a causal inverse, but if the 
> > entire signal can be buffered before processing,
> > then it seems to me that the inverse does not
> > need to be causal. Am I missing something?
> > 
> > In article <f56893ae.0408190125.27ec4942@posting.google.com>,
> > Rune Allnor <allnor@tele.ntnu.no> wrote:
> > >shafik@u.arizona.edu (Shafik) wrote in message news:<5e5cc96b.0408182053.6c5bab12@posting.google.com>...
> > >> Hello everyone,
> > >> 
> > >>    Can anyone tell me if there is a known way to find the inverse of a
> > >> given FIR filter? The idea is that given some signal y(x), and an FIR
> > >> filter f(x), I want to
> > >> find a way to compute g(x) which is the inverse of f(x) such that:
> > >> 
> > >>       y(x) * f(x) * g(x) = y(x).  
> > >>       (* = convolution)
> > >> 
> > >> Meaning if I have some signal, I can apply an FIR filter to shape it
> > >> then apply an inverse FIR to bring back the original signal.
> > >> 
> > >> Any ideas?
> > >
> > >You can get what you want, provided the FIR filter is minimum phase, i.e.
> > >has all its zeros in the interior of the unit circle. In that case, you 
> > >can get a stable, causal IIR filter that compensates the effects of 
> > >your FIR filter. 
> > >
> > >See e.g. section 4.6 in
> > >
> > >Proakis & Manolakis: "Digital Signal processing. Pronciples, Algorithms
> > >   and Applications" Prentice-Hall, 3rd ed. 1996.
> > >
> > >Rune

Shafik wrote:

> Hello everyone,
> 
>    Can anyone tell me if there is a known way to find the inverse of a
> given FIR filter?

Generalized Levinson-Durbin in Matlab:


function inv=invimplms(den,n,d)
%syntax	inv=invimplms(den,n,d)
%	den - denominator FIR (to be inverted)
%	n   - length of result
%	d   - delay of result
%	inv - inverse FIR of length n with delay d
%
% Author: Bob Cain, May 1, 2001 arcane@arcanemethods.com

	m=xcorr(den,n-1);
	m=m(n:end);
	b=[den(d+1:-1:1);zeros(n-d-1,1)];
	inv=toepsolve(m,b);

%----------------------------------------------------

function quo=divimplms(num,den,n,d)
%syntax	quo=divimplms(num,den,n,d)
%	num - numerator FIR
%	den - denominator FIR
%	n   - length of result
%	d   - delay of result
%	quo - quotient FIR of length n delayed by d
%
% Author: Bob Cain, May 1, 2001 arcane@arcanemethods.com

	m=xcorr(den,n-1);
	m=m(n:end);
	b=xcorr([zeros(d,1);num],den,n-1);
	b=b(n:-1:1);
	quo=toepsolve(m,b);

%----------------------------------------------------

function hinv=toepsolve(r,q)
% Solve Toeplitz system of equations
%	R*hinv = q
%		where R is the symmetric Toeplitz matrix
%		who's first column is r
%
%	Assumes all inputs are real
%
%	Inputs:
%		r - first column of Toeplitz matrix, length n
%		q - rhs vector, length n
%	Outputs:
%		hinv - length n solution
%
%   Algorithm from Roberts & Mullis, p.233
%
%   Author: T. Krauss, Sept 10, 1997

	n=size(q);
	a=zeros(n+1,2);
	a(1,1)=1;c=1;
	hinv=zeros(n+1,1);
	hinv(1)=q(1)/r(1);
	alpha=r(1);
	d=2;

	for k = 1:n,
	   a(k+1,c)=0;
	   a(1,d)=1;
	   beta=0;
	   j=1:k;
	   beta=sum(r(k+2-j).*a(j,c)')/alpha;
	   a(j+1,d)=a(j+1,c)-beta*a(k+1-j,c);
	   alpha = alpha*(1 - beta^2);
	   hinv(k+1,1) = ( q(k+1)-sum(r(k+2-j).*hinv(j,1)') ) / alpha;
	   hinv(j) = hinv(j)+a(k+2-j,d)*hinv(k+1);
	   temp=c;
	   c=d;
	   d=temp;
	end

%----------------------------------------------------

And toepsolve in C for compilation in Matlab:

#include <math.h>
#include "mex.h"

/*  function h = toepsolve(r,q);
  *  TOEPSOLVE  Solve Toeplitz system of equations.
  *    Solves R*h = q, where R is the symmetric Toeplitz matrix
  *    whos first column is r
  *    Assumes all inputs are real
  *    Inputs:
  *       r - first column of Toeplitz matrix, length n
  *       q - rhs vector, length n
  *    Outputs:
  *       h - length n solution
  *
  *   Algorithm from Roberts & Mullis, p.233
  *
  *   Author: T. Krauss, Sept 10, 1997
  */
void mexFunction(
     int nlhs,
     mxArray *plhs[],
     int nrhs,
     const mxArray *prhs[]
)
{
    double *a,*h,beta;
    int j,k;

    double eps = mxGetEps();
    int n = (mxGetN(prhs[0])>=mxGetM(prhs[0])) ? mxGetN(prhs[0]) : mxGetM(prhs[0]) ;
    double *r = mxGetPr(prhs[0]);
    double *q = mxGetPr(prhs[1]);
    double alpha = r[0];

    n = n - 1;

    plhs[0] = mxCreateDoubleMatrix(n+1,1,0);
    h = mxGetPr(plhs[0]);

    h[0] = q[0]/r[0];

    a = mxCalloc((n+1)*(n+1),sizeof(double));
    if (a == NULL) {
        mexErrMsgTxt("Sorry, failed to allocate buffer.");
    }

    a[(0*(n+1))+0] = 1.0;

    for (k = 1; k <= n; k++) {
        a[(k*(n+1))+k-1] = 0;
        a[(0*(n+1))+k] = 1.0;
        beta = 0.0;
        for (j = 0; j <= k-1; j++) {
            beta += r[k-j]*a[(j*(n+1))+k-1];
        }
        beta /= alpha;
        for (j = 1; j <= k; j++) {
            a[(j*(n+1))+k] = a[(j*(n+1))+k-1] - beta*a[((k-j)*(n+1))+k-1];
        }
        alpha *= (1 - beta*beta);
        h[k] = q[k];
        for (j = 0; j <= k-1; j++) {
            h[k] -= r[k-j]*h[j];
        }
        h[k] /= alpha;
        for (j = 0; j <= k-1; j++) {
            h[j] += a[((k-j)*(n+1))+k]*h[k];
        }
    } /* loop over k */

    mxFree(a);

    return;
}


toepsolve() is a memory hog and it doesn't need to be.
It blows a circulant matrix up to full size when modular
indexing on one row could be done.  Never got around to
fixing it and I should because it severely limits the
size problem it will handle.


Bob
-- 

"Things should be described as simply as possible, but no simpler."

                                              A. Einstein