fircls1.m

有关matlab的电子书籍有一定的帮助希望有用

function [h,a] = fircls1(n,wo,dp,ds,varargin)
%FIRCLS1 Low & high pass FIR filter design by constrained least-squares.
%   B = FIRCLS1(N,WO,DP,DS) is a length N+1 linear-phase lowpass FIR filter
%   with cut-off frequency W0 and maximum band deviations DP and DS. 
%   W0 is between 0 and 1 (1 corresponds to half the sampling frequency),
%   DP is the maximum passband deviation from 1, and DS is the maximum stopband
%   deviation from 0.
% 
%   B = FIRCLS1(N,WO,DP,DS,'high') is a highpass filter (order N must be even).
%
%   B = FIRCLS1(N,WO,DP,DS,WT) and B = FIRCLS1(N,WO,DP,DS,WT,'high') meets
%   a passband or stopband edge requirement. If WT is in the passband, then 
%   the use of this parameter ensures that |E(WT)| <= DP where E(w) is the 
%   error function. Similarly, if WT is in the stopband, then the use of WT
%   ensures that |E(WT)| <= DS. Note that in the design of very narrow band 
%   filters with small DP and DS, there may not exist a filter of the given 
%   length that meets these specifications. 
%
%   B = FIRCLS1(N,WO,DP,DS,WP,WS,K) weights the square error in the passband 
%   K times greater than that in the stopband. WP is the passband edge of the 
%   L2 weight function and WS is the stopband edge (WP < W0 < WS).  Use 
%   trailing WT and 'high' arguments to meet a requirement or design a highpass
%   filter as in the case with no weighting function, e.g. 
%      FIRCLS1(N,WO,DP,DS,WP,WS,K,WT,'high'). 
%   In the 'high' pass filter case, you must have WS < W0 < WP.
%
%   EXAMPLES
%      n = 55;
%      wo = 0.3;
%      dp = 0.02; ds = 0.008;
%      h = fircls1(n,wo,dp,ds);           % no L2 weights
%      wp = 0.28; ws = 0.32;
%      K = 10;
%      h1 = fircls1(n,wo,dp,ds,wp,ws,K);  % L2 weight
%
%   MONITORING THE DESIGN
%      For a textual progress report on the iteration, use a trailing 'trace' 
%      argument, e.g. FIRCLS1(N,W0,DP,DS,...,'trace').  To display plots of 
%      the design iteration, use a trailing 'plots' argument, 
%      FIRCLS1(N,W0,DP,DS,...,'plots').  For both text and plots, use 'both'.
%
%   See also FIRCLS, 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
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.11 $  $Date: 1997/12/02 18:37:01 $

%  subprograms : local_max.m, frefine.m

% --------- check input paramters ---------

LOW = 0;
HIGH = 1;
pass_type = LOW;
EDGE = 0;
WEIGHTS = 0;
TEXT_PF = 0;
PLOT_PF = 0;

if nargin >= 7 & ~isstr(varargin{3})
	WEIGHTS = 1;
	wp = varargin{1};
	ws = varargin{2};
	K = varargin{3};
	varargin = varargin(4:end);
end

if nargin < 4
	error('You did not supply enough input paramters.')
elseif ( wo <= 0) | ( wo >= 1)
	error('W0 must lie between 0 and 1.')
elseif prod([size(n), size(wo), size(dp), size(ds)]) > 1
	error('Each of N, WO, DP, and DS must be a scalar.')
elseif (dp <= 0) | (ds <= 0)
	error('Both DP and DS must be positive.')
elseif length(varargin)==1
	if isstr(varargin{1})
		switch lower(varargin{1}(1))
		case 'h'
			pass_type = HIGH;
		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 string input.')
		end
	else
		wt = varargin{1};
		EDGE = 1;
	end
elseif length(varargin)==2
	if isstr(varargin{1})
               	switch lower(varargin{1}(1))
               	case 'h'
                       	pass_type = HIGH;
               	otherwise
                       	error('Unrecognized string input.')
               	end
               	switch lower(varargin{2}(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 string input.')
               	end
	else
		wt = varargin{1};
		EDGE = 1;
               	switch lower(varargin{2}(1))
               	case 'h'
                       	pass_type = HIGH;
		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 string input.')
               	end
	end
elseif length(varargin)==3
	wt = varargin{1};
	EDGE = 1;
	if isstr(varargin{2})
                switch lower(varargin{2}(1))
                case 'h'
                        pass_type = HIGH;
                otherwise
                        error('Unrecognized string input.')
                end
	end
	if isstr(varargin{3})
                switch lower(varargin{3}(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 string input.')
                end
	end
end
PASS = 1;
STOP = 2;
if EDGE
   if prod(size(wt)) > 1
	error('WT must be a scalar')
   elseif (wt <= 0) | (wt >= 1) 
	error('WT must lie between 0 and 1.')
   elseif wt == wo
	error('WT must not equal W0')
   elseif wt < wo
	if pass_type == LOW
		edge_type = PASS;
	else
		edge_type = STOP;
	end
   else
	if pass_type == LOW
		edge_type = STOP;
	else
		edge_type = PASS;
	end
   end
end

if WEIGHTS
    if pass_type == LOW
            if (wp > wo) | (ws < wo)
                    error('For lowpass filters, WP < W0 < WS is needed.')
            end
    else
            if (ws > wo) | (wp < wo)
                    error('For highpass filters, WS < W0 < WP is needed.')
            end
    end
end

if pass_type == LOW
   up = [1+dp  ds]; d1 = dp;
   lo = [1-dp -ds]; d2 = ds;
   mag = [1 0];
else
   up = [ ds 1+dp]; d1 = ds;
   lo = [-ds 1-dp]; d2 = dp;
   mag = [0 1];
end

n = n+1;  % convert order to length for this routine
if rem(n,2) == 1
	parity = 1;
else
	parity = 0;
	if pass_type == HIGH
		error('High pass filters must have EVEN orders.')
	end
end

% --------- start algorithm ---------

wo = wo*pi;
if WEIGHTS
    wp = wp*pi;
    ws = ws*pi;
end
if EDGE, wt = wt*pi; end
L = 2^ceil(log2(5*n));
r = sqrt(2);
w = [0:L]*pi/L;                  % w includes both 0 and pi
q = round(wo*L/pi);
u = [up(1)*ones(q,1); up(2)*ones(L+1-q,1)];
l = [lo(1)*ones(q,1); lo(2)*ones(L+1-q,1)];
if parity  == 1
	m = (n-1)/2;
	if WEIGHTS
		c = 2*[wp/r; [sin(wp*[1:m])./[1:m]]']/pi;
	else
		c = 2*[wo/r; [sin(wo*[1:m])./[1:m]]']/pi;
	end
	if pass_type == HIGH
		c = -c;
		c(1) = c(1)+r;
	end
	Z = zeros(2*L-n,1);
	v = [0:m];
	tt = 1 - 1/r;
	NP = m+1;	% NP : number of parameters
else
	m = n/2;
        v = [1:m]-1/2;
	% c = [4*((sin(wo*[2*[1:m]-1]/2))./(2*[1:m]-1))/pi]';
	if WEIGHTS
		c = [2*sin(wp*v)./(v*pi)]';
	else
        	c = [2*sin(wo*v)./(v*pi)]';
	end
	Z = zeros(4*L,1);
	tt = 0;
	NP = m;		% NP : number of parameters
end


if WEIGHTS
	% ----- construct R matrix --------------------
	if rem(n,2) == 1
   	% odd length symmetric filter
   		v1 = 1:m;
   		v2 = m:2*m;
   		if pass_type == LOW
      			tp = [wp+K*(pi-ws), (sin(v1*wp)-K*sin(v1*ws))./v1]/pi;
      			hk = ((sin(v2*wp)-K*sin(v2*ws))./v2)/pi;
   		else % pass_type == HIGH
      			tp = [(pi-wp)+K*ws, (-sin(v1*wp)+K*sin(v1*ws))./v1]/pi;
      			hk = ((-sin(v2*wp)+K*sin(v2*ws))./v2)/pi;
   		end
   		R = toeplitz(tp,tp) + hankel(tp,hk);
   		R(1,:) = R(1,:)/r;
   		R(:,1) = R(:,1)/r;
   		Ri = inv(R);
	else
   	% even length symmetric filter
   		v1 = 1:(m-1);
   		v2 = (m-1):2*(m-1);
   		tp = [(wp+K*(pi-ws)), (sin(v1*wp)-K*sin(v1*ws))./v1]/pi;
   		v1 = 1:m;
   		v2 = m:2*m-1;
   		tp2 = ((sin(v1*wp)-K*sin(v1*ws))./v1)/pi;
   		hk2 = ((sin(v2*wp)-K*sin(v2*ws))./v2)/pi;
   		R = toeplitz(tp,tp) + hankel(tp2,hk2);
   		Ri = inv(R);
	end

	a = Ri*c;       % best L2 cosine coefficients
else
	a = c;          % best L2 cosine coefficients
end
mu = [];                % Lagrange multipliers
SN = 1e-8;              % Small Number
% -------- BEGIN LOOPING --------------
while 1

   % calculate H 
   if parity == 1
      H = fft([a(1)*r; a(2:m+1); Z; a(m+1:-1:2)]); 
      H = real(H(1:L+1)/2);
   else
      Z(2:2:2*m) = a;
      Z(4*L-2*m+2:2:4*L) = a(m:-1:1);
      H = fft(Z);
      H = real(H(1:L+1)/2);
   end

   % find local maxima and minima
   kmax = local_max(H);
   kmin = local_max(-H);
   % if filter length is even, then remove w=pi from constraint set
   if parity == 0
        n1 = length(kmax);
        if kmax(n1) == L+1, kmax(n1) = []; n1 = n1 - 1; end
        n2 = length(kmin);
        if kmin(n2) == L+1, kmin(n2) = []; n2 = n2 - 1; end
   end

   % refine frequencies
   cmax = (kmax-1)*pi/L;      cmin = (kmin-1)*pi/L;
   cmax = frefine(a,v,cmax);  cmin = frefine(a,v,cmin);

   % insert wt into constraint set if neccessary
   if EDGE
   if pass_type == LOW
      if edge_type == PASS
         w_key = max(cmax(cmax<wo));
         if wt > w_key
            kmin = [kmin; 1];
            cmin = [cmin; wt];
         end
      else % edge_type == STOP
         w_key = min(cmin(cmin>wo));
         if wt < w_key
            kmax = [kmax; L];
            cmax = [cmax; wt];
         end
      end
   else % pass_type == HIGH
      if edge_type == PASS
         w_key = min(cmax(cmax>wo));
         if wt < w_key
            kmin = [kmin; L];
            cmin = [cmin; wt];
         end
      else % edge_type == STOP
         w_key = max(cmin(cmin<wo));
         if wt > w_key
            kmax = [kmax; 1];
            cmax = [cmax; wt];
         end
      end
   end
   end

   % evaluate H at refined frequencies
   Hmax = cos(cmax*v)*a - tt*a(1);
   Hmin = cos(cmin*v)*a - tt*a(1);

   % determine new constraint set
   v1   = Hmax > u(kmax)-100*SN;
   v2   = Hmin < l(kmin)+100*SN;
   kmax = kmax(v1); kmin = kmin(v2);
   cmax = cmax(v1); cmin = cmin(v2);
   Hmax = Hmax(v1); Hmin = Hmin(v2);
   n1   = length(kmax);
   n2   = length(kmin);

   % plot figures
   if PLOT_PF
   wv = [cmax; cmin];
   Hv = [Hmax; Hmin];
	subplot(311)
        plot(w/pi,H),
	axis([0 1 -.2 1.2])
	hold on,
	plot(wv/pi,Hv,'o'),
	hold off
	subplot(312)
	plot(w/pi,H-mag(1)),
	hold on,
	plot(wv/pi,Hv-mag(1),'o'),
	hold off
	axis([0 wo/pi -2*d1 2*d1])
	subplot(313)
	plot(w/pi,H-mag(2)),
	hold on,
	plot(wv/pi,Hv-mag(2),'o'),
	hold off
	axis([wo/pi 1 -2*d2 2*d2])
   pause(.5)
   end

   % remove a constraint set frequency if neccessary (if otherwise overdetermined)
   if (n1+n2) > NP
        if parity == 1
           H0 =  a(1)/sqrt(2) + sum(a(2:m+1));
           Hpi = a(1)/sqrt(2) + sum(a(3:2:m+1)) - sum(a(2:2:m+1));
           dH0dw = -sum(a(2:m+1)'.*([1:m].^2));
           dHpidw = sum(a(2:2:m+1)'.*([1:2:m].^2)) - sum(a(3:2:m+1)'.*([2:2:m].^2));
           if dH0dw > 0, E0 = lo(1) - H0;
           else, E0 = H0 - up(1); end
           if dHpidw > 0, Epi = lo(2) - Hpi;
           else, Epi = Hpi - up(2); end
        else % parity == 0
		% when length is even, we know that
		%  we must remove w = 0;
		E0 = 0; Epi = 1;
        end
        if E0 > Epi
                % remove w = pi from constraint set
                [temp1, k1] = max(kmin);
                [temp2, k2] = max(kmax);
                if temp1 < temp2
                        % pi is in kmax
                        kmax(k2) = [];
                        cmax(k2) = [];
                        Hmax(k2) = [];
                        n1 = n1 - 1;
                else
                        % pi is in kmin
                        kmin(k1) = [];
                        cmin(k1) = [];
                        Hmin(k1) = [];
                        n2 = n2 - 1;
                end
        else
		% remove w = 0 from constraint set
                [temp1, k1] = min(kmin);
                [temp2, k2] = min(kmax);
                if temp1 < temp2
                        % 0 is in kmin
                        kmin(k1) = []; 
			cmin(k1) = [];
                        Hmin(k1) = [];
                        n2 = n2 - 1;
                else
                        % 0 is in kmax
                        kmax(k2) = [];
                        cmax(k2) = [];
                        Hmax(k2) = [];
                        n1 = n1 - 1;
                end
        end
   end

   % check stopping criterion
   E  = max([Hmax-u(kmax); l(kmin)-Hmin; 0]);
   if TEXT_PF
      fprintf(1,'    Bound Violation = %15.13f  \n',E);
   end
   if E < SN
      break
   end

   % calculate new Lagrange multipliers
   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]*([1:m]-1/2));
   end
   d = [u(kmax); -l(kmin)];

   if WEIGHTS
       mu = (G*Ri*G')\(G*Ri*c-d);
   else
       mu = (G*G')\(G*c-d);            
   end

   % iteratively remove negative multiplier
   [min_mu,K] = min(mu);
   while min_mu < 0
      G(K,:) = [];
      d(K) = [];
      if WEIGHTS
          mu = (G*Ri*G')\(G*Ri*c-d);
      else
          mu = (G*G')\(G*c-d);            
      end
      [min_mu,K] = min(mu);
   end

   % determine new cosine coefficients
   if WEIGHTS
      a = Ri*(c-G'*mu);
   else
      a = c-G'*mu;
   end

end

if parity == 1
   h = [a(m+1:-1:2)/2; a(1)/r; a(2:m+1)/2]';
else
   h = [a(m:-1:1); a]'/2;
end

if nargout > 1
    a = 1;
end