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