Complex Least Squares FIR filter design using Levinson' s algorithm

Question

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;
end
a=a/L;
b=b/L;
% 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));
end

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

om1=pi*[linspace(0,0.23,230),linspace(0.3,0.5,200),linspace(0.57,1,430)];
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

Jeff

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

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.