i'm trying to implement an arbitrary IIR filter design. I'm decomposing all the algorithm into smaller pieces that i can test.

I currently am on the part where it designs a FIR filter from an arbitrary response.

Matlab's main function is here :

function h = lslevin (N ,om, D, W)
% Complex Least Squares FIR filter design using Levinson' s algorithm
% h filter impulse response
% N filter length
% om frequency grid (0 <= om <= pi)
% D complex desired frequenccy response on the grid om
% W positive weighting function on the grid om
om = om(:);
D = D(:); 
W = W(:); 
L = length(om);
DR = real(D) ;
DI = imag(D);
 a = zeros(N,1); 
 b = a;
% Set up vectors for quadratic objective function
% (avoid building matrices)
dvec = D;
evec = ones(L, 1); 
e1 = exp(j*om);
for i=1:N,
    a(i) = W.'*real (evec);
    b(i) = W.'*real (dvec);
    evec = evec.*e1; 
    dvec = dvec.*e1;
% Compute weighted l2 solution
h = levin(a, b) ;

and levin function :

function x = levin(a ,b)
% function x = levin(a ,b)
% solves system of complex linear equations toeplitz (a)*x = b
% using Levinson's algorithm
% a ... first row of positive definite Hermitian Toeplitz matrix
% b ... right hand side vector
% Author: Mathias C. Lang, Vienna University of Technology, AUSTRIA
% 9-97
a = a(:);
b = b(:); 
n = length(a);
t = 1; 
alpha = a(1); 
x = b(1) / a(1);
for i = 1:n-1,
    k = -(a(i+1 :-1: 2)' * t)/alpha;
    t = [t;0] + k*flipud( [conj(t);0]);
    alpha = alpha*( 1 - abs(k)^2);
    k = (b(i+1) - a(i+1:-1: 2)' * x) /alpha;
    x = [x;0] + k*flipud(conj(t));

an example function can be done here (it should be a simple bandpass):

D1=[zeros(1,230), exp(-j*om1(231:430) *2), zeros(1,430)]; % zeros are bandstop others bandpass
W1=[10* ones(1,230), ones(1,200), 10 *ones(1,430) ];%weight function
h = lslevin(61,om1, D1, W1) ;
wavwrite(h, 44100,  'e:\tmp2.wav'); % write the impulse response

but when i compute the Impulse Response, i have something like that :

which is far from the original design

can you help me pinpoint the errors in the code ?

thanks a lot


EDIT1 : here is the final FFT of the IR from the answer given below:


1 Answer 1


From this term in D1:

exp(-j*om1(231:430) *2)

it looks like you're asking it for a filter with a group delay of 2 samples even though the number of taps is 61. Try setting this to:

exp(-j*om1(231:430) * 61/2)

i.e. set a desired group delay of about half the filter length.

  • 1
    $\begingroup$ oooh it works! amazed $\endgroup$
    – ion_one
    Aug 28, 2013 at 14:22
  • $\begingroup$ Woo-hoo! Glad to be of help. Any chance of appending the results to the question? It might help someone in the future. $\endgroup$
    – Peter K.
    Aug 28, 2013 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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