| << Chapter < Page | Chapter >> Page > |
The following is the author's implementation of Soewito's linearization. This program is based on the theoretical development presented in [link] , and is designed to be used by other programs in this work. For the original implementation (which includes a number of additions to handle numerical issues) the reader should consult the original work by A. Soewito [link] .
function [bb,aa] = l2soe(b,a,D,varargin)% This program is an implementation of Soewito's quasilinearization
% to design least squares filters using a solution error criterion.% This implementation accepts arbitrary samples between 0 and pi,
% and requires an initial filter guess.%
% [B,A]= L2SOE(Bo,Ao,D) designs an optimal L2 approximation to a
% complex response D that is specified over a uniform frequency% grid between 0 and PI, using as initial point the vector {Bo,Ao}.
%% [B,A] = L2SOE(Bo,Ao,D,F) designs an optimal L2 approximation to a% complex response D that is specified over a nonuniform frequency
% grid F defined between 0 and PI.%
% [B,A]= L2SOE(Bo,Ao,D,F,W) designs a weighted L2 approximation,
% with weighting function W.%
% [B,A]= L2SOE(Bo,Ao,D,F,W,ITER) finds a solution in at most ITER
% iterations (the default for this value is 30 iterations).%
% [B,A]= L2SOE(Bo,Ao,D,F,W,ITER,TOL) defines convergence when the
% norm of the updating vector is less than TOL (default is 0.01).%
% [B,A]= L2SOE(Bo,Ao,D,...,'diags') plots the desired and current
% spectra at each iteration and displays error measurements.% Author: R. A. Vargas, 4-20-08error(nargchk(3,8,nargin))if isempty(varargin), fdiags = false;
elsefdiags = varargin(length(varargin));
if (ischar(fdiags{1})&&strcmp(fdiags,'diags')),
fdiags = true; varargin(length(varargin)) = [];
else fdiags = false; endend
if length(varargin)<4 % pad varargin with []'svarargin{4} = [];end
[f,W,mxit,etol]= deal(varargin{:});
if isempty(W), W = ones(size(D)); end; W = W(:);if isempty(mxit), mxit = 30; end
if isempty(etol), etol = 0.01; endM = length(b); N = length(a); % Numerator and denominator lengthsd = D(:); % Form conjugate symmetric desired response and weights
if (~(isreal(d(1))&&isreal(d(length(d))) )),
error('Real filters have real spectra values at 0 and Pi');end
D = [d; flipud(conj(d(2:length(d)-1)))];
W = [W; flipud(conj(W(2:length(W)-1)))];
% Define frequencies over whole unit circleif isempty(f), f = [0:2*pi/length(D):2*pi*(length(D)-1)/length(D)]';else f = [f; 2*pi-flipud(f(2:length(f)-1))]; end% Initialize relevant variables
h = [b; a(2:N)];
F = [exp(-i.*(f*[0:M-1])) -diag(D)*exp(-i.*(f*[1:N-1]))];
b0 = zeros(size(b)); a0 = zeros(size(a)); c = 0;% Iterative loopwhile (norm([b0;a0]-[b;a],2)>= etol&&c<mxit),
c = c+1; b0 = b; a0 = a;% Vector update
V = diag(W.*freqz(1,a,f));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))];
% Diagnosticsif fdiags,
sprintf(strcat('Iteration = %g,',' Error norm = %g '),...c,norm(D-freqz(b,a,f),2))
[hh,ff]= freqz(b,a,1024,'whole');
plot(f./(2*pi),abs(D),'.',ff./(2*pi),abs(hh))title(strcat('Iteration:',int2str(c))); pause
endend%if fdiags,
if c==mxit, disp('Maximum number of L2 iterations reached')else sprintf('L2 convergence reached after %g iterations',c)
endbb = real(b); aa = real(a);
Notification Switch
Would you like to follow the 'Iterative design of l_p digital filters' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|