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
EDIT1 : here is the final FFT of the IR from the answer given below: