This program implements the algorithm presented in to design an [Math Processing Error] magnitude IIR filter. The program combines ideas used in the programs mentioned above, including [Math Processing Error]-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[Math Processing Error] algorithms to generate a suitable initial guess for the magnitude [Math Processing Error] 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 parameters
M = 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 p
K = 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 indexes
ih = [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 Prony
mxit = 100; % Max. iterations for Prony stage
etol = 0.01; % Error tolerance for Prony stage
k = 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 update
cM = 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 p
while (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 p
cP = 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 update
cM = 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 Update
F = [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 Quasilinearization
while (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 Update
b = (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
Posts
0 comments:
Post a Comment