Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# mattdsp/levin.m

Last active Jan 2, 2018
Least squares FIR digital filter design with prescribed magnitude AND phase response using Levinson's algorithm
 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 % 1997-09 % mattsdspblog@gmail.com 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
 function h = lslevin(N,om,D,W) % 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 frequency response on the grid om % W positive weighting function on the grid om % % example: length 61 bandpass, band edges [.23,.3,.5,.57]*pi, % weighting 1 in passband and 10 in stopbands, desired passband % group delay 20 samples % % om=pi*[linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)]; % D=[zeros(1,230),exp(-j*om(231:430)*20),zeros(1,430)]; % W=[10*ones(1,230),ones(1,200),10*ones(1,430)]; % h = lslevin(61,om,D,W); % % Author: Mathias C. Lang, Vienna University of Technology % 1998-07 % mattsdspblog@gmail.com 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);
to join this conversation on GitHub. Already have an account? Sign in to comment