This program designs a constrained [Math Processing Error] 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
x = find(fd>ft); Ad(x) = zeros(size(x)); % Add zeros on stopband
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);
e = C*a - Ad;
end
%%% Recovery of h from a
a(1) = 2*a(1);
h = [flipud(a(2:length(a)));a]./2;
Posts
0 comments:
Post a Comment