| << Chapter < Page | Chapter >> Page > |
This section includes Matlab implementations of some of the algorithms described in this work.
The following program designs a Type-I length-
linear phase
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
[link] .
% Lp Linear Phase design program. Uses p-homotopy and partial updating.% Author: R. A. Vargas, 4-20-08%%% User parameters
P = 400; % Desired pK = 1.7; % Homotopy rate
fp = 0.2; fs = 0.24; % Normalized passband&stopband
L = 21; % Filter length%%% InitializationNF = 2^ceil(log2(10*L)); % Number of frequency samples
Np = round(NF*fp/(.5 - (fs - fp))); % No. of passband samplesdp = fp/(Np - 1);
Ns = NF - Np; % No. of stopband samplesds = (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 guesse = 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 updateW = diag(w/sum(w)); % Normalized weights
x = (W*C)\(W*Ad); % WLS solutiona = (x + (p-2)*a) ./(p-1); % Partial update
e = C*a - Ad; % Errorend%%% Recovery of h from a
a(1) = 2*a(1);h = [flipud(a(2:length(a)));a]./2;
The following code expands on the program from
"The linear phase
FIR algorithm" . If at a given iteration the error increases, the idea is to take a step back in
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 parametersP = 200; % Desired p
K = 1.7; % Homotopy ratefp = 0.2; fs = 0.248; % Normalized passband&stopband
L = 21; % Filter lengthdk = 0.9; % K update factor%%% Initialization
NF = 2^ceil(log2(10*L)); % Number of frequency samplesNp = 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]'; % FrequenciesAd = [ones(Np,1); zeros(Ns,1)]; % Desired responseM = floor((L+1)/2)-1;
C = cos(2*pi*fd*(0:M)); % Fourier matrixa = C\Ad; a0 = a; % L2 initial guess
e = C*a - Ad; en = e; % Initial error%%% Algorithm iterationc = 1; maxiter = 40; p = 2;
while (c<maxiter),
p = min(P,K*p); % p-homotopy updatew = abs(e).^((p-2)/2); % Lp weight update
W = diag(w/sum(w)); % Normalized weightsx = (W*C)\(W*Ad); % WLS solution
a = (x + (p-2)*a) ./(p-1); % Partial updateen = 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 stepend
end%%% Recovery of h from aa(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;
Notification Switch
Would you like to follow the 'Iterative design of l_p digital filters' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|