The following program designs an [Math Processing Error] IIR filter given a complex-valued desired frequency response. Note the similarities of the program compared to the FIR ones.
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 lengths
w = 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 frequencies
D = [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 quasilinearization
b = (bb + (p-2)*b) ./(p-1); % Partial update for b
a = (aa + (p-2)*a) ./(p-1); % Partial update for a
end
Posts
0 comments:
Post a Comment