gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\fircls.m
function [h,a] = fircls(n,be,mag,up,lo,PF)
%FIRCLS Linear-phase FIR filter design by constrained least-squares.
% B = FIRCLS(N,F,A,UP,LO) is a length N+1 linear phase FIR filter which
% approximates the desired piecewise constant frequency response in F and A
% subject to the constraints given in UP and LO.
%
% A is a vector describing the piecewise constant desired amplitude of
% the frequency response. The length of A is the number of different bands.
% UP and LO are vectors of the same length of A. They give the upper and
% lower bounds for the frequency response in each band.
% F is a vector of transition frequencies. These frequencies must
% be in increasing order, must start with 0.0, and must end with 1.0.
% The length of F should exceed the length of A by exactly 1.
%
% EXAMPLE
% n = 51;
% f = [0 0.4 0.8 1];
% a = [0 1 0];
% up = [ 0.02 1.02 0.01];
% lo = [-0.02 0.98 -0.01];
% b = fircls(n,f,a,up,lo);
%
% NOTE
% By setting LO equal to 0 in the stopbands, a nonnegative frequency
% response amplitude can be obtained. Such filters can be spectrally
% factored to obtain minimum phase filters.
%
% MONITORING THE DESIGN
% For a textual progress report on the iteration, use a trailing 'trace'
% argument, e.g. FIRCLS(N,F,A,UP,LO,'trace'). To display plots of
% the design iteration, use a trailing 'plots' argument,
% FIRCLS(N,F,A,UP,LO,'plots'). For both text and plots, use 'both'.
%
% See also FIRCLS1, FIRLS, REMEZ.
% Reference:
% I. W. Selesnick, M. Lang and C. S. Burrus,
% Constrained Least Square Design of FIR Filters
% Without Specified Transition Bands,
% IEEE Trans. on Signal Processing,
% Aug 1996, Vol 44, Number 8, pp. 1879-1892.
% Author: Ivan Selesnick, Rice University, 1995
% Copyright (c) 1988-98 by The MathWorks, Inc.
% $Revision: 1.1 $ $Date: 1998/06/03 14:42:39 $
% subprograms : local_max.m, frefine.m
% -------- check input parameters -------------
if nargin < 5
error('You did not supply enough input parameters.')
else
lm = length(mag);
lu = length(up);
ll = length(lo);
if (lm ~= lu) | (lm ~= ll) | (lu ~= ll)
error('Each of A, UP, and LO must have the same length.')
elseif length(be) ~= (lm+1)
error('F must have one more element than A')
elseif be(1) ~= 0
error('F must begin with 0')
elseif be(lm+1) ~= 1
error('F must end with 1')
elseif any(diff(be)<=0)
error('The frequencies in F must be in ascending order.')
elseif any(up <= lo)
error('Each element of UP must be greater than the corresponding element of LO.')
elseif any(lo > mag)
error('Each entry of LO should no greater than the corresponding element of A.')
elseif any(up < mag)
error('Each entry of UP should no less than the corresponding element of A.')
end
end
n = n+1; % convert order to length
TEXT_PF = 0;
PLOT_PF = 0;
if nargin == 6
PF = lower(PF);
switch PF(1)
case 't' % as in 'text'
TEXT_PF = 1;
case 'p' % as in 'plots'
PLOT_PF = 1;
case 'b' % as in 'both'
TEXT_PF = 1; PLOT_PF = 1;
otherwise,
error('Unrecognized trace option. Try ''trace'', ''plots'', or ''both''.')
end
end
if prod(size(n)) > 1
error('n should be a scalar')
end
L = 2^ceil(log2(3*n));
num_bands = length(mag);
if rem(n,2) == 0
parity = 0;
m = n/2;
if (lo(num_bands) > 0) | (up(num_bands) < 0)
error('Even length filters must equal 0 at F=1.');
end
else
parity = 1;
m = (n-1)/2;
end
r = sqrt(2);
be = be * pi;
% ----- calculate Fourier coefficients and upper ---
% ----- and lower bound functions ------------------
if parity == 1
c = zeros(m+1,1);
Z = zeros(2*L-1-2*m,1);
v = [0:m];
tt = 1-1/r;
else
c = zeros(m,1);
Z = zeros(4*L,1);
v = [1:m]-1/2;
tt = 0;
end
u = [];
l = [];
for k = 1:num_bands
if parity == 1
c = c + mag(k) * ...
(2*[be(k+1)/r; [sin(be(k+1)*[1:m])./[1:m]]']/pi - ...
2*[be(k)/r; [sin(be(k)*[1:m])./[1:m]]']/pi);
else
c = c + mag(k) * ...
( [4*((sin(be(k+1)*[2*[1:m]-1]/2))./(2*[1:m]-1))/pi]' - ...
[4*((sin(be(k) * [2*[1:m]-1]/2))./(2*[1:m]-1))/pi]' ) ;
end
q = round(L*(be(k+1)-be(k))/pi);
u = [u; up(k)*ones(q,1)];
l = [l; lo(k)*ones(q,1)];
end
flen = length(u);
if flen < L+1
ov = ones(L+1-flen,1);
u(flen+1:L+1) = up(num_bands)*ov;
l(flen+1:L+1) = lo(num_bands)*ov;
elseif flen > L+1
u = u(1:L+1);
l = l(1:L+1);
end
w = [0:L]'*pi/L;
a = c; % best L2 cosine coefficients
mu = []; % Lagrange multipliers
SN = 1e-8; % Small Number
it = 0;
kmax = zeros(0,1); kmin = zeros(0,1);
cmax = zeros(0,1); cmin = zeros(0,1);
okmax = []; okmin = []; uvo = 0;
ocmax = []; ocmin = []; lvo = 0;
while 1
if (uvo < -SN/10) | (lvo < -SN/10)
% ----- include old extremal ----------------
if uvo < lvo
kmax = [kmax; okmax(k1)]; okmax(k1) = [];
cmax = [cmax; ocmax(k1)]; ocmax(k1) = [];
else
kmin = [kmin; okmin(k2)]; okmin(k2) = [];
cmin = [cmin; ocmin(k2)]; ocmin(k2) = [];
end
else
% ----- calculate A -------------------------
if parity == 1
A = fft([a(1)*r;a(2:m+1);Z;a(m+1:-1:2)]);
A = real(A(1:L+1))/2;
else
Z(2:2:2*m) = a;
Z(4*L-2*m+2:2:4*L) = a(m:-1:1);
A = fft(Z);
A = real(A(1:L+1)/2);
end
if PLOT_PF
cc = [cmax; cmin];
occ = [ocmax; ocmin];
if ~isempty(cc)
Acc = cos(cc*v)*a - tt*a(1);
else
Acc = [];
end
if ~isempty(occ)
Aocc = cos(occ*v)*a - tt*a(1);
else
Aocc = [];
end
subplot(length(be),1,1)
plot(w/pi,A), hold on
plot(cc/pi,Acc,'o')
plot(occ/pi,Aocc,'x')
hold off
for k = 1:length(be)-1
subplot(length(be),1,k+1);
plot(w/pi,A), hold on
plot(cc/pi,Acc,'o')
plot(occ/pi,Aocc,'x')
hold off
axis([be(k)/pi be(k+1)/pi ...
mag(k)-1.5*(mag(k)-lo(k)) mag(k)+1.5*(up(k)-mag(k))])
ylabel(sprintf('Band #%g',k))
end
xlabel('Frequency')
drawnow
end
% ----- find extremals ----------------------
okmax = kmax; okmin = kmin;
ocmax = cmax; ocmin = cmin;
kmax = local_max(A); kmin = local_max(-A);
if parity == 0
n1 = length(kmax);
if kmax(n1) == L+1, kmax(n1) = []; end
n2 = length(kmin);
if kmin(n2) == L+1, kmin(n2) = []; end
end
cmax = (kmax-1)*pi/L; cmin = (kmin-1)*pi/L;
cmax = frefine(a,v,cmax); cmin = frefine(a,v,cmin);
Amax = cos(cmax*v)*a - tt*a(1);
Amin = cos(cmin*v)*a - tt*a(1);
v1 = Amax > u(kmax)-SN;
v2 = Amin < l(kmin)+SN;
kmax = kmax(v1); kmin = kmin(v2);
cmax = cmax(v1); cmin = cmin(v2);
Amax = Amax(v1); Amin = Amin(v2);
% ----- check stopping criterion ------------
Eup = Amax-u(kmax); Elo = l(kmin)-Amin;
E = max([Eup; Elo]);
if TEXT_PF
fprintf(1,' Bound Violation = %15.13f \n',E);
end
if E < SN, break, end
end
% ----- calculate new multipliers -----------
n1 = length(kmax); n2 = length(kmin);
if parity == 1
G = [ones(n1,m+1); -ones(n2,m+1)] .* cos([cmax;cmin]*[0:m]);
G(:,1) = G(:,1)/r;
else
G = [ones(n1,m); -ones(n2,m)] .* cos([cmax;cmin]*v);
end
d = [u(kmax); -l(kmin)];
mu = (G*G')\(G*c-d);
% ----- remove negative multiplier ----------
[min_mu,K] = min(mu);
while min_mu < 0
G(K,:) = []; d(K) = [];
mu = (G*G')\(G*c-d);
if K > n1
kmin(K-n1) = [];
cmin(K-n1) = [];
%if length(Amin)>=K-n1
% Amin(K-n1) = [];
%end
n2 = n2 - 1;
else
kmax(K) = [];
cmax(K) = [];
%Amax(K) = [];
n1 = n1 - 1;
end
[min_mu,K] = min(mu);
end
% ----- determine new coefficients ----------
a = c-G'*mu;
if length(ocmax)>0
oAmax = cos(ocmax*v)*a - tt*a(1);
[uvo,k1] = min(u(okmax) - oAmax);
else
uvo = 0;
end
if length(ocmin)>0
oAmin = cos(ocmin*v)*a - tt*a(1);
[lvo,k2] = min(oAmin - l(okmin));
else
lvo = 0;
end
it = it + 1;
end
if parity == 1
h = [a(m+1:-1:2); a(1)*r; a(2:m+1)]/2;
else
h = [a(m:-1:1); a]/2;
end
h = h(:)';
if nargout > 1
a = 1;
end