implementation of Soewito's quasilinearization

The following is the author's implementation of Soewito's linearization. This program is based on the theoretical development presented in [1], 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 [1].

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-08
 
error(nargchk(3,8,nargin))
if isempty(varargin), fdiags = false;
else
   fdiags = varargin(length(varargin));
   if (ischar(fdiags{1}) && strcmp(fdiags,'diags')),
      fdiags = true; varargin(length(varargin)) = [];
   else fdiags = false; end
end
if length(varargin)<4    % pad varargin with []'s
   varargin{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; end
 
M = length(b); N = length(a); % Numerator and denominator lengths
d = 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 circle
if 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 loop
while (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))];
   % Diagnostics
   if 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
   end
end
 
%if fdiags,
if c==mxit, disp('Maximum number of L2 iterations reached')
else sprintf('L2 convergence reached after %g iterations',c)
end
bb = real(b); aa = real(a);