# The magnitude [Math Processing Error] IIR algorithm

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 parametersM = 5; N = 5;                % Numerator & denominator lengthsfp = 0.2; fs = 0.24;         % Normalized passband & stopbandt = 129;                     % Number of samples between 0 & piP = 30;                      % Desired pK = 1.4;                     % Homotopy rate %%% Initializationw = [0:pi/(t-1):pi]';        % Radial frequencyip = max(find(w<=fp*2*pi));  % Passband indexesis = min(find(w>=fs*2*pi));  % Stopband indexesit = 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 fh(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 stageetol = 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 stagemxitM = 50; mxit2 = 50;% Error tolerances for Solution Error L2 stageetolM = 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))];   endend %%% Magnitude Lp Iterative Method% Max. iterations for Solution Error Lp stagemxitP = 200; mxitM = 60; mxit2 = 50;% Error tolerances for Solution Error Lp stageetolP = 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 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``

This post first appeared on Matlab Projects....., please read the originial post: here

# Share the post

The magnitude [Math Processing Error] IIR algorithm

×