remezord.m

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

function [N, ff, aa, wts] = remezord(fcuts, mags, devs, fsamp, cellflag)
%REMEZORD  FIR order estimator (lowpass, highpass, bandpass, multiband)
%   [N,Fo,Ao,W] = REMEZORD(F,A,DEV,Fs) finds the approximate order N, 
%   normalized frequency band edges Fo, frequency band magnitudes Ao and 
%   weights W to be used by the REMEZ function as follows:
%       B = REMEZ(N,Fo,Ao,W)
%   The resulting filter will approximately meet the specifications given
%   by the input parameters F, A, and DEV.  F is a vector of cutoff
%   frequencies in Hz, in ascending order between 0 and half the sampling
%   frequency Fs. If you do not specify Fs, it defaults to 2.  A is a vector
%   specifying the desired function's amplitude on the bands defined by F.
%   The length of F is twice the length of A, minus 2 (it must therefore
%   be even).  The first frequency band always starts at zero, and the last
%   always ends at Fs/2.  It is not necessary to add these elements to the 
%   F vector.  DEV is a vector of maximum deviations or ripples allowable 
%   for each band.
%
%   C = REMEZORD(F,A,DEV,FS,'cell') is a cell-array whose elements are the 
%   parameters to REMEZ.
%
%   EXAMPLE
%      Design a lowpass filter with a passband cutoff of 1500, a 
%      stopband cutoff of 2000Hz, passband ripple of 0.01, stopband ripple 
%      of 0.1, and a sampling frequency of 8000Hz:
%
%         [n,fo,mo,w] = remezord( [1500 2000], [1 0], [0.01 0.1], 8000 );
%         b = remez(n,fo,mo,w);
%
%      This is equivalent to
%         c = remezord( [1500 2000], [1 0], [0.01 0.1], 8000, 'cell');
%         b = remez(c{:});
%
%   CAUTION 1: The order N is often underestimated. If the filter does not
%   meet the original specifications, a higher order such as N+1 or N+2 will.
%   CAUTION 2: Results are inaccurate if cutoff frequencies are near zero
%   frequency or the Nyquist frequency.
%
%   See also REMEZ, KAISERORD.

%   Author(s): J. H. McClellan, 10-28-91
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.15 $  $Date: 1997/11/26 20:13:30 $

%   References:
%     [1] Rabiner & Gold, Theory and Applications of DSP, pp. 156-7.     

error(nargchk(3,5,nargin))
typ = 'Bandpass';
if nargin == 3,
    fsamp = 2;
end
fcuts = fcuts/fsamp;       %  NORMALIZE to sampling frequency

% Turn vectors into column vectors
fcuts = fcuts(:);
mags = mags(:);
devs = devs(:);

[mf,nf] = size(fcuts);
[mm,nm] = size(mags);
[md,nd] = size(devs);
nbands = mm;

if length(mags) ~= length(devs)
    error('Requires M and DEV to be the same length.')
end

if mf ~= 2*(nbands-1)
    error('Length of F must be 2*length(M)-2.')
end

zz = mags==0;             % find stopbands
devs = devs./(zz+mags);   % divide delta by mag to get relative deviation

% Determine the smallest width transition zone
% Separate the passband and stopband edges
%
f1 = fcuts(1:2:(mf-1));
f2 = fcuts(2:2:mf);
[df,n] = min(f2-f1);

%=== LOWPASS case: Use formula (ref: Herrmann, Rabiner, Chan)
%
if( nbands==2 )
     L = remlpord( f1(n), f2(n), devs(1), devs(2));

%=== BANDPASS case:
%    - try different lowpasses and take the WORST one that
%        goes thru the BP specs; try both transition widths
%    - will also do the bandreject case
%    - does the multi-band case, one bandpass at a time.
%    
else
  L = 0;
  for i=2:nbands-1,
    L1 = remlpord( f1(i-1), f2(i-1), devs(i),   devs(i-1) );
    L2 = remlpord( f1(i),   f2(i),   devs(i),   devs(i+1) );
    L = max( [L; L1; L2] );
  end
end

N = ceil( L ) - 1;   % need order, not length, for remez

%=== Make the MATLAB compatible specs for remez
%
ff = [0;2*fcuts;1];
am(1:2:2*nbands-1) = mags;	
aa = [am';0] + [0;am'];
wts = ones(size(devs))*max(devs)./devs;

if nargout == 1 & nargin == 5
   N = {N, ff, aa, wts};
end

function L = remlpord(freq1, freq2, delta1, delta2 )
%REMLPORD FIR lowpass filter Length estimator
%
%   L = remlpord(freq1, freq2, dev1, dev2)
%
%   input:
%     freq1: passband cutoff freq (NORMALIZED)
%     freq2: stopband cutoff freq (NORMALIZED)
%     dev1:  passband ripple (DESIRED)
%     dev2:  stopband attenuation (not in dB)
%
%   output:
%     L = filter Length (# of samples)
%         NOT the order N, which is N = L-1
%
%   NOTE: Will also work for highpass filters (i.e., f1 > f2)
%         Will not work well if transition zone is near f = 0, or
%         near f = fs/2

% 	Author(s): J. H. McClellan, 10-28-91
%       Was Revision: 1.4, Date: 1994/01/25 17:59:46
	
%   References:
%     [1] Rabiner & Gold, Theory and Appications of DSP, pp. 156-7.     

AA = [-4.278e-01  -4.761e-01  0;
      -5.941e-01   7.114e-02  0;
      -2.660e-03   5.309e-03  0 ];
d1 = log10(delta1);
d2 = log10(delta2);
D = [1 d1 d1*d1] * AA * [1; d2; d2*d2];
%
bb = [11.01217; 0.51244];
fK =  [1.0  (d1-d2)] * bb;
%
df = abs(freq2 - freq1);
%
L = D/df - fK*df + 1;