| << Chapter < Page | Chapter >> Page > |
This program designs a constrained linear phase filter without requiring a fixed transition band. Based on a transition frequency the program determines at each iteration an induced transition band and creates weights to "constrain" the error. One of the key steps is determining the frequencies that exceed the constrains and, from these frequencies, determine the induced transition band.
% Code for CLS design. This example designs a linear phase filter.
% No transition bands are described.% Author: R. A. Vargas, 4-20-08%%% User parametersP = 80; % Desired p
K = 1.7; % p-homotopy rate
ft = 0.25; % Transition frequencyL = 21; % Filter length
tol = 0.02; % Constraint tolerance%%% InitializationNF = 2^ceil(log2(10*L)); % No. of frequency samples
fd = linspace(0,.5,NF)'; % Linearly spaced sampled frequenciesAd = ones(NF,1); % Desired frequency response
x = find(fd>ft); Ad(x) = zeros(size(x)); % Add zeros on stopband
C = cos(2*pi*fd*[0:floor((L+1)/2)-1]); % Fourier matrix
a = C\Ad; % L2 initial guesse = C*a - Ad; % Initial error%%% Algorithm iteration
i = 60; % Maximum number of iterationsc = 1; p = 2;
for m = 1:i, m,pp = min(K*p,P); % P-homotopy
w = 1 + abs((e./(0.95*tol)).^((p-2)/2)); % Polynomial weightingX = local_max(-abs(e)); % Find all local extrema
% Here we find the index of the two extrema near the trans bandEp = max(X(find(fd(X)<ft))); % Passband extrema
Es = min(X(find(fd(X)>ft))); % Stopband extrema
w(Ep:Es) = 1; % Unweighting of trans band% WLS solution w/partial update
a = ((p-2)*a + ((diag(w)*C)\(diag(w)*Ad))) ./ (p-1);e = C*a - Ad;
end%%% Recovery of h from aa(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;
The following program designs an
IIR filter given a complex-valued desired frequency response. Note the similarities of the program compared to the FIR ones. This program calls the function
l2soe from
"An implementation of Soewito's quasilinearization" to solve the weighted nonlinear
problem.
function [b,a] = cplx_lp_iir(b0,a0,D,w,P);% function [b,a] = CPLX_LP_IIR(b0,a0,D,w,p);% This function designs an Lp IIR filter given a complex desired
% frequency response arbitrarily sampled between 0 and PI. The% algorithm used is an Iterative Reweighted Least Squares (IRLS)
% method. For the WLS part, a quasilinearization L2 IIR algorithm% is used.
%% [b0,a0] : Initial guess% D : Desired complex response (defined between 0 and PI)
% w : Frequency samples between 0 and PI% p : Desired maximum p% Author: R. A. Vargas, 4-20-08%%% Initialization
[b,a]= l2soe(b0,a0,D,w); % Initial guess
M = length(b); N = length(a); % Numerator and denominator lengthsw = w(:); d = D(:);
if (~(isreal(d(1))&&isreal(d(length(d))) )),
error('Real filters have real spectra values at 0 and Pi');end
% Form conjugate symmetric desired response and frequenciesD = [d; flipud(conj(d(2:length(d)-1)))];f = [w; 2*pi-flipud(w(2:length(w)-1))];K = 1.7; p = 2;
mxit = 100; etol = 0.01; c = 0;b0 = zeros(size(b)); a0 = zeros(size(a));% Algorithm iteration
while ((norm([b0;a0]-[b;a],2)>= etol&c<mxit) | p<P),
c = c+1; b0 = b; a0 = a;p = min(P,K*p); % p homotopy
W = abs(freqz(b,a,w) - d).^((p-2)/2);[bb,aa] = l2soe(b,a,d,w,W,60); % L2 quasilinearizationb = (bb + (p-2)*b) ./(p-1); % Partial update for b
a = (aa + (p-2)*a) ./(p-1); % Partial update for aend
Notification Switch
Would you like to follow the 'Iterative design of l_p digital filters' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|