# Connexions

You are here: Home » Content » Iterative Design of l_p Digital Filters » Appendix: Matlab Code

### Recently Viewed

This feature requires Javascript to be enabled.

Inside Collection:

Collection by: Ricardo Vargas. E-mail the author

# Appendix: Matlab Code

Module by: Ricardo Vargas. E-mail the author

This section includes Matlab implementations of some of the algorithms described in this work.

## The linear phase lplp FIR algorithm

The following program designs a Type-I length-LL linear phase lplp filter with fixed transition bands. The code creates even sampling in the bands but it can easily be modified to accomodate for uneven sampling. For other linear phase types the user can modify the definition of C according to (Reference).

% Lp Linear Phase design program. Uses p-homotopy and partial updating.

% Author: R. A. Vargas, 4-20-08

%%% User parameters
P = 400;                         % Desired p
K = 1.7;                         % Homotopy rate
fp = 0.2; fs = 0.24;             % Normalized passband & stopband
L  = 21;                         % Filter length

%%% Initialization
NF = 2^ceil(log2(10*L));         % Number of frequency samples
Np = round(NF*fp/(.5 - (fs - fp))); % No. of passband samples
dp = fp/(Np - 1);
Ns = NF - Np;                    % No. of stopband samples
ds = (0.5 - fs)/(Ns - 1);
fd = [(0:Np - 2)*dp,fp,(0:Ns - 2)*ds + fs,0.5]'; % Frequencies
Ad = [ones(Np,1); zeros(Ns,1)];  % Desired response
M = floor((L+1)/2)-1;
C = cos(2*pi*fd*(0:M));          % Fourier matrix
a = C\Ad;                        % L2 initial guess
e = C*a - Ad;                    % Initial error

%%% Algorithm iteration
c = 1; i = 40; p = 2;
while (c<i),
c = c+1;
p = min(P,K*p);              % p-homotopy update
w = abs(e).^((p-2)/2);       % Lp weight update
W = diag(w/sum(w));          % Normalized weights
x = (W*C)\(W*Ad);            % WLS solution
a = (x + (p-2)*a) ./(p-1);   % Partial update
e = C*a - Ad;                % Error
end

%%% Recovery of h from a
a(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;


## The adaptive linear phase lplp FIR algorithm

The following code expands on the program from "The linear phase lplp FIR algorithm". If at a given iteration the error increases, the idea is to take a step back in pp and then take a smaller step. This is achieved by reducing the homotopy step parameter K in the program.

% Lp Linear Phase design program. Uses p-homotopy and partial updating.
% This code uses the adaptive algorithm.

% Author: R. A. Vargas, 4-20-08

%%% User parameters
P = 200;                         % Desired p
K = 1.7;                         % Homotopy rate
fp = 0.2; fs = 0.248;            % Normalized passband & stopband
L  = 21;                         % Filter length
dk = 0.9;                        % K update factor

%%% Initialization
NF = 2^ceil(log2(10*L));         % Number of frequency samples
Np = round(NF*fp/(.5 - (fs - fp))); % No. of passband samples
dp = fp/(Np - 1);
Ns = NF - Np;                    % No. of stopband samples
ds = (0.5 - fs)/(Ns - 1);
fd = [(0:Np - 2)*dp,fp,(0:Ns - 2)*ds + fs,0.5]'; % Frequencies
Ad = [ones(Np,1); zeros(Ns,1)];  % Desired response
M = floor((L+1)/2)-1;
C = cos(2*pi*fd*(0:M));          % Fourier matrix
a = C\Ad; a0 = a;                % L2 initial guess
e = C*a - Ad; en = e;            % Initial error

%%% Algorithm iteration
c = 1; maxiter = 40; p = 2;
while (c<maxiter),
p = min(P,K*p);               % p-homotopy update
w = abs(e).^((p-2)/2);        % Lp weight update
W = diag(w/sum(w));           % Normalized weights
x = (W*C)\(W*Ad);             % WLS solution
a = (x + (p-2)*a) ./(p-1);    % Partial update
en = C*a - Ad;                % Error
if (norm(en,P) <= norm(e,P)) | p>=P
c = c+1;
e = en;
a0 = a;
else
a = a0;
p = p/K;
K = dk*K;                  % Update homotopy step
end
end

%%% Recovery of h from a
a(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;


## The constrained l2l2 FIR algorithm

This program designs a constrained l2l2 linear phase filter without requiring a fixed transition band. Based on a transition frequency the program determines at each iteration an induced transition band and creates weights to "constrain" the error. One of the key steps is determining the frequencies that exceed the constrains and, from these frequencies, determine the induced transition band.

% Code for CLS design. This example designs a linear phase filter.
% No transition bands are described.

% Author: R. A. Vargas, 4-20-08

%%% User parameters
P = 80;                    % Desired p
K = 1.7;                   % p-homotopy rate
ft = 0.25;                 % Transition frequency
L = 21;                    % Filter length
tol = 0.02;                % Constraint tolerance

%%% Initialization
NF = 2^ceil(log2(10*L));   % No. of frequency samples
fd = linspace(0,.5,NF)';   % Linearly spaced sampled frequencies
Ad = ones(NF,1);           % Desired frequency response
C = cos(2*pi*fd*[0:floor((L+1)/2)-1]);    % Fourier matrix
a = C\Ad;                  % L2 initial guess
e = C*a - Ad;              % Initial error

%%% Algorithm iteration
i = 60;                    % Maximum number of iterations
c = 1; p = 2;
for m = 1:i, m,p
p = min(K*p,P);                     % P-homotopy
w = 1 + abs((e./(0.95*tol)).^((p-2)/2));  % Polynomial weighting
X = local_max(-abs(e));             % Find all local extrema
% Here we find the index of the two extrema near the trans band
Ep = max(X(find(fd(X)<ft)));        % Passband extrema
Es = min(X(find(fd(X)>ft)));        % Stopband extrema
w(Ep:Es) = 1;                       % Unweighting of trans band
% WLS solution w/partial update
a = ((p-2)*a + ((diag(w)*C)\(diag(w)*Ad))) ./ (p-1);
end

%%% Recovery of h from a
a(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;


## The complex lplp IIR algorithm

The following program designs an lplp IIR filter given a complex-valued desired frequency response. Note the similarities of the program compared to the FIR ones. This program calls the function l2soe from "An implementation of Soewito's quasilinearization" to solve the weighted nonlinear l2l2 problem.

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


## An 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);


## The magnitude lplp IIR algorithm

This program implements the algorithm presented in (Reference) to design an lplp magnitude IIR filter. The program combines ideas used in the programs mentioned above, including pp-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 l2l2 algorithms to generate a suitable initial guess for the magnitude lplp 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


## References

1. Soewito, Atmadji W. (1990, December). Least Square Digital Filter Design in the Frequency Domain. Ph. D. Thesis. Rice University.

## Content actions

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

#### Collection to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

#### Module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks