gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\remezord.m
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.1 $ $Date: 1998/06/03 14:43:40 $ % 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;