| << Chapter < Page | Chapter >> Page > |
This program implements the algorithm presented in
[link] to design an
magnitude IIR filter. The program combines ideas used in the programs mentioned above, including
-homotopy and partial filter updating, together with the concept of phase updating (to achieve magnitude approximation) and the quasilinearization implementation from
"An implementation of Soewito's quasilinearization" . To improve on the convergence properties, the program initially implements equation error and solution error
algorithms to generate a suitable initial guess for the magnitude
iteration. it is important torealize that each of the three computatino al stages of this code can be fine-tuned in terms of convergence parameters. This program uses even sampling but it can easily be modified for uneven sampling (by merely changing the initial definition of
f ). the program also defines a ramp transition band but the user can define any desired real response by defining
h below.
% This program designs a magnitude lp IIR digital filter.
% Data is uniformly sampled between 0 and PI.% Author: R. A. Vargas, 4-20-08%%% User parametersM = 5; N = 5; % Numerator&denominator lengths
fp = 0.2; fs = 0.24; % Normalized passband&stopband
t = 129; % Number of samples between 0&pi
P = 30; % Desired pK = 1.4; % Homotopy rate%%% Initialization
w = [0:pi/(t-1):pi]'; % Radial frequency
ip = max(find(w<=fp*2*pi)); % Passband indexes
is = min(find(w>=fs*2*pi)); % Stopband indexes
it = ip+1:is-1; % Transition band indexesih = [1:ip is:t-1]; % Indexes at which error is measured% Form conjugate symmetric desired response D and frequency f
h(1:ip) = ones(ip,1); h(is:t) = zeros(t-is+1,1);h(ip+1:is-1) = ((w(ip+1:is-1)/2/pi)-fs)./(fp-fs); d = h(:);
D = [d; flipud(conj(d(2:length(d)-1)))];
f = [w; 2*pi-flipud(w(2:length(w)-1))];
L = length(D); % Number of samples on unit circle%%% Equation Error Magnitude L2 estimation via Pronymxit = 100; % Max. iterations for Prony stage
etol = 0.01; % Error tolerance for Prony stagek = 2*etol; c = 0; a0 = zeros(N,1); b0 = zeros(M,1);
while (k>=etol&&c<mxit),
c = c+1;h = ifft(D);
H = h(toeplitz([1:L]',[1 L:-1:L-N+2]));H1 = H(1:M,:); h2 = H(M+1:L,1);
H2 = H(M+1:L,2:size(H,2));a = [1; -H2\h2];b = H1*a;
k = norm([b0;a0]-[b;a],2);a0 = a; b0 = b;
G = fft(b,L)./fft(a,L);D = abs(D).*G./abs(G);
end%%% Solution Error Magnitude L2 estimation via Quasilinearization% Max. iterations for Solution Error L2 stage
mxitM = 50; mxit2 = 50;% Error tolerances for Solution Error L2 stage
etolM = 0.005; etol2 = 0.01;cM = 0; bM = zeros(size(b)); aM = zeros(size(a));
while (norm([bM;aM]-[b;a],2)>=etolM&&cM<mxitM),
% Initialize relevant variables at each phase updatecM = cM+1; bM = b; aM = a;
G = fft(b,L)./fft(a,L);D = abs(D).*G./abs(G); % Phase Update
h = [b; a(2:N)];
F = [exp(-i.*(f*[0:M-1])) -diag(D)*exp(-i.*(f*[1:N-1]))];
b2 = zeros(size(b)); a2 = zeros(size(a)); c2 = 0;%%% Complex L2 loop using Quasilinearization
while (norm([b2;a2]-[b;a],2)>=etol2&&c2<mxit2),
c2 = c2+1; b2 = b; a2 = a;V = diag(freqz(1,a,f)); % Vector update
H = freqz(b,a,f);G = [exp(-i.*(f*[0:M-1])) -diag(H)*exp(-i.*(f*[1:N-1]))];h = (V*G)\(V*(D+H-F*h));
b = h(1:M);a = [1; h(M+1:length(h))];end
end%%% Magnitude Lp Iterative Method% Max. iterations for Solution Error Lp stage
mxitP = 200; mxitM = 60; mxit2 = 50;% Error tolerances for Solution Error Lp stage
etolP = 0.005; etolM = 0.005; etol2 = 0.005;W = ones(size(d)); W = [W; flipud(conj(W(2:length(W)-1)))];bP = zeros(size(b)); aP = zeros(size(a)); cP = 0; p = 2*K;
%%% Outer loop to update pwhile (norm([bP;aP]-[b;a],2)>= etolP&&cP<mxitP&&p<=P),
if p>=P, etolP = 0.0001; end
% Initialize relevant variables at each update of pcP = cP + 1; bP = b; aP = a;
bM = zeros(size(b)); aM = zeros(size(a)); cM = 0;%%% Magnitude Lp loop via phase update
while (norm([bM;aM]-[b;a],2)>= etolM&&cM<mxitM),
% Initialize relevant variables at each phase updatecM = cM+1; bM = b; aM = a;
h = [b; a(2:N)];
b2 = zeros(size(b)); a2 = zeros(size(a)); c2 = 0;G = freqz(b,a,f);
D = abs(D).*G./abs(G); % Phase UpdateF = [exp(-i.*(f*[0:M-1])) -diag(D)*exp(-i.*(f*[1:N-1]))];E = abs(D - freqz(b,a,f));
W = E.^((p-2)/2);W(it) = W(it)./4;
%%% Complex Lp loop via WCL2 using Quasilinearizationwhile (norm([b2;a2]-[b;a],2)>= etol2&&c2<mxit2),
c2 = c2+1; b2 = b; a2 = a;V = diag(W.*freqz(1,a,f)); % Vector update
H = freqz(b,a,f);G = [exp(-i.*(f*[0:M-1])) -diag(H)*exp(-i.*(f*[1:N-1]))];h = (V*G)\(V*(D+H-F*h));
b = h(1:M);a = [1; h(M+1:length(h))];end
% Partial Updateb = (b + (p-2)*bM) ./(p-1);
a = (a + (p-2)*aM) ./(p-1);end
G = fft(b,L)./fft(a,L);D = abs(D).*G./abs(G);
p = min(P,K*p);end
Notification Switch
Would you like to follow the 'Iterative design of l_p digital filters' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|