|
回复 2楼 的帖子
我的matlabR14(V7.04)里确实有fir1和fir2,在Signal Processing Toolbox里。eight院长没有可能是matlab的版本问题吧。另外,楼主,关于二者的区别就在help文档的Note里写的很清楚了啊。
Note Use fir2 for windowed filters with arbitrary frequency response.
Note Use fir1 for windows-based standard lowpass, bandpass, highpass, and bandstop configurations
fir1
Wwindow-based finite impulse response filter design
Syntax
b = fir1(n,Wn)
b = fir1(n,Wn,'ftype')
b = fir1(n,Wn,window)
b = fir1(n,Wn,'ftype',window)
b = fir1(...,'normalization')
Description
fir1 implements the classical method of windowed linear-phase FIR digital filter design [1]. It designs filters in standard lowpass, highpass, bandpass, and bandstop configurations. By default the filter is normalized so that the magnitude response of the filter at the center frequency of the passband is 0 dB.
Note Use fir2 for windowed filters with arbitrary frequency response. fir2
Ffrequency sampling-based finite impulse response filter design
Syntax
b = fir2(n,f,m)
b = fir2(n,f,m,window)
b = fir2(n,f,m,npt)
b = fir2(n,f,m,npt,window)
b = fir2(n,f,m,npt,lap)
b = fir2(n,f,m,npt,lap,window)
Description
fir2 designs frequency sampling-based digital FIR filters with arbitrarily shaped frequency response.
Note Use fir1 for windows-based standard lowpass, bandpass, highpass, and bandstop configurations.
fir1.m- function [b,a] = fir1(N,Wn,varargin)
- %FIR1 FIR filter design using the window method.
- % B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter
- % and returns the filter coefficients in length N+1 vector B.
- % The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0
- % corresponding to half the sample rate. The filter B is real and
- % has linear phase. The normalized gain of the filter at Wn is
- % -6 dB.
- %
- % B = FIR1(N,Wn,'high') designs an N'th order highpass filter.
- % You can also use B = FIR1(N,Wn,'low') to design a lowpass filter.
- %
- % If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an
- % order N bandpass filter with passband W1 < W < W2. You can
- % also specify B = FIR1(N,Wn,'bandpass'). If Wn = [W1 W2],
- % B = FIR1(N,Wn,'stop') will design a bandstop filter.
- %
- % If Wn is a multi-element vector,
- % Wn = [W1 W2 W3 W4 W5 ... WN],
- % FIR1 returns an order N multiband filter with bands
- % 0 < W < W1, W1 < W < W2, ..., WN < W < 1.
- % B = FIR1(N,Wn,'DC-1') makes the first band a passband.
- % B = FIR1(N,Wn,'DC-0') makes the first band a stopband.
- %
- % B = FIR1(N,Wn,WIN) designs an N-th order FIR filter using
- % the N+1 length vector WIN to window the impulse response.
- % If empty or omitted, FIR1 uses a Hamming window of length N+1.
- % For a complete list of available windows, see the help for the
- % WINDOW function. KAISER and CHEBWIN can be specified with an
- % optional trailing argument. For example, B = FIR1(N,Wn,kaiser(N+1,4))
- % uses a Kaiser window with beta=4. B = FIR1(N,Wn,'high',chebwin(N+1,R))
- % uses a Chebyshev window with R decibels of relative sidelobe
- % attenuation.
- %
- % For filters with a gain other than zero at Fs/2, e.g., highpass
- % and bandstop filters, N must be even. Otherwise, N will be
- % incremented by one. In this case the window length should be
- % specified as N+2.
- %
- % By default, the filter is scaled so the center of the first pass band
- % has magnitude exactly one after windowing. Use a trailing 'noscale'
- % argument to prevent this scaling, e.g. B = FIR1(N,Wn,'noscale'),
- % B = FIR1(N,Wn,'high','noscale'), B = FIR1(N,Wn,wind,'noscale'). You
- % can also specify the scaling explicitly, e.g. FIR1(N,Wn,'scale'), etc.
- %
- % See also KAISERORD, FIRCLS1, FIR2, FIRLS, FIRCLS, CFIRPM,
- % FIRPM, FREQZ, FILTER, WINDOW.
- % FIR1 is an M-file implementation of program 5.2 in the IEEE
- % Programs for Digital Signal Processing tape.
- % Author(s): L. Shure
- % L. Shure, 4-5-90, revised
- % T. Krauss, 3-5-96, revised
- % Copyright 1988-2004 The MathWorks, Inc.
- % $Revision: 1.15.4.3 $ $Date: 2004/10/18 21:08:06 $
- % Reference(s):
- % [1] "Programs for Digital Signal Processing", IEEE Press
- % John Wiley & Sons, 1979, pg. 5.2-1.
- error(nargchk(2,5,nargin));
- % Parse optional input arguments
- [Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin{:});
- error(msg);
- % Compute the frequency vector
- [nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype);
- error(msg);
- % Compute the magnitude vector
- [aa,First_Band] = desiredmag(Ftype,nbands);
- % Check for appropriate filter order, increase when necessary
- [N,msg1,msg2] = firchk(N,ff(end),aa);
- error(msg1);
- warning(msg2);
- % Work with filter length (= order + 1)
- L = N + 1;
- % Check for valid window, or assign default if empty
- [Wind,msg] = chkwindow(Wind,L);
- error(msg);
- % Compute unwindowed impulse response
- hh = firls(L-1,ff,aa);
- % Window impulse response to get the filter
- b = hh.*Wind(:)';
- a = 1;
- if SCALING,
- % Scale so that passband is approx 1
- b = scalefilter(b,First_Band,ff,L);
- end
- %-------------------------------------------------------------------------
- function [Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin)
- %PARSEOPTARGS Parse optional input arguments.
- % Up to 3 optional input arguments, always in this order:
- % 1 - Filter type flag, can be 'low','high','bandpass','stop','DC-0','DC-1'
- % 2 - Window vector
- % 3 - 'noscale' flag
- % Initialize output args.
- SCALING = [];
- [Ftype,Wind,Scale] = assignoptargs(Wn,varargin{:});
- [Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale);
- if ~isempty(msg),
- return
- end
- switch lower(Scale),
- case 'noscale';
- SCALING = 0;
- case 'scale';
- SCALING = 1;
- end
- %--------------------------------------------------------------------------
- function [Ftype,Wind,Scale] = assignoptargs(Wn,varargin)
- %ASSIGNOPTARGS Assign optional input arguments to the appropriate variables.
- % default optional parameter values:
- Wind = [];
- Scale = 'scale';
- Ftype = defaultftype(Wn);
- switch length(varargin)
- case 1
- if ischar(varargin{1}) && (length(varargin{1})>0),
- s = upper(varargin{1});
- switch upper(s)
- case {'SCALE','NOSCALE'}
- Scale = s;
- otherwise
- Ftype = s;
- end
- else
- Wind = varargin{1};
- end
- case 2
- if ischar(varargin{1})
- Ftype = varargin{1};
- else
- Wind = varargin{1};
- end
- if ischar(varargin{2})
- Scale = varargin{2};
- else
- Wind = varargin{2};
- end
- case 3
- Ftype = varargin{1};
- Wind = varargin{2};
- Scale = varargin{3};
- end
- %--------------------------------------------------------------------------
- function [Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale)
- %VALIDATEARGS Test if arguments are valid.
- msg = '';
- % Assign a default Ftype when an empty is given. Backwards compatibility
- if isempty(Ftype),
- Ftype = defaultftype(Wn);
- end
- Ftypeopts = {'LOW','HIGH','BANDPASS','STOP','DC-0','DC-1'};
- Scaleopts = {'NOSCALE','SCALE'};
- indx = find(strncmpi(Ftype, Ftypeopts, length(Ftype)));
- if isempty(indx),
- msg = 'Unrecognized or ambiguous filter type specified.';
- return
- else
- Ftype = Ftypeopts{indx};
- end
- scaleindx = find(strncmpi(Scale, Scaleopts, length(Scale)));
- if isempty(scaleindx),
- msg = 'Scaling option must be ''noscale'' or ''scale''.';
- return
- else
- Scale = Scaleopts{scaleindx};
- end
- if ~any(size(Wind) <= 1),
- msg = 'The window specified must be a vector';
- return
- else
- Wind = Wind(:).'; % Make it a row vector
- end
- %--------------------------------------------------------------------------
- function [nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype)
- %DESIREDFREQ Compute the vector of frequencies to pass to FIRLS.
- %
- % Inputs:
- % Wn - vector of cutoff frequencies.
- % Ftype - string with desired response ('low','high',...)
- %
- % Outputs:
- % nbands - number of frequency bands.
- % ff - vector of frequencies to pass to FIRLS.
- % Ftype - converted filter type (if it's necessary to convert)
- % Initialize output args.
- nbands = [];
- ff = [];
- msg = '';
- if any( Wn<0 | Wn>1 ),
- msg = 'Frequencies must fall in range between 0 and 1.';
- return
- end
- if any(diff(Wn)<0),
- msg = 'Frequencies must be increasing';
- return
- end
- Wn = Wn(:)';
- nbands = length(Wn) + 1;
- if (nbands > 2) && strcmpi(Ftype,'bandpass'),
- Ftype = 'DC-0'; % make sure default 3 band filter is bandpass
- end
- ff = [0,Wn(1:nbands-1); Wn(1:nbands-1),1];
- ff = ff(:);
- %-------------------------------------------------------------------------
- function [aa,First_Band] = desiredmag(Ftype,nbands)
- %DESIREDMAG Compute the magnitude vector to pass to FIRLS.
- First_Band = isempty(findstr('DC-0',Ftype)) && isempty(findstr('HIGH',Ftype));
- mags = rem( First_Band + (0:nbands-1), 2);
- aa = [mags(:)'; mags(:)'];
- aa = aa(:);
- %--------------------------------------------------------------------------
- function [Wind,msg] = chkwindow(Wind,L)
- %CHKWINDOW Check if specified window is valid, assign default if empty.
- msg = '';
- if isempty(Wind),
- % Replace the following with the default window of your choice.
- Wind = hamming(L);
- end
-
- if length(Wind) ~= L
- msg = 'The window length must be the same as the filter length.';
- end
- %
- % to use Kaiser window, beta must be supplied
- % att = 60; % dB of attenuation desired in sidelobe
- % beta = 0.1102*(att-8.7);
- % wind = kaiser(L,beta);
- %---------------------------------------------------------------------------
- function b = scalefilter(b,First_Band,ff,L)
- %SCALEFILTER Scale fitler to have passband approx. equal to one.
- if First_Band
- b = b / sum(b); % unity gain at DC
- else
- if ff(4)==1
- % unity gain at Fs/2
- f0 = 1;
- else
- % unity gain at center of first passband
- f0 = mean(ff(3:4));
- end
- b = b / abs( exp(-j*2*pi*(0:L-1)*(f0/2))*(b.') );
- end
- %----------------------------------------------------------------------------
- function Ftype = defaultftype(Wn)
- %DEFAULTFTYPE Assign default filter type depending on number of bands.
- if length(Wn) == 1,
- Ftype = 'low';
- elseif length(Wn) == 2,
- Ftype = 'bandpass';
- elseif length(Wn) >= 3,
- Ftype = 'dc-0';
- end
复制代码 fir2.m- function [b,a] = fir2(nn, ff, aa, npt, lap, wind)
- %FIR2 FIR arbitrary shape filter design using the frequency sampling method.
- % B = FIR2(N,F,A) designs an N'th order FIR digital filter with the
- % frequency response specified by vectors F and A, and returns the
- % filter coefficients in length N+1 vector B. Vectors F and A specify
- % the frequency and magnitude breakpoints for the filter such that
- % PLOT(F,A) would show a plot of the desired frequency response.
- % The frequencies in F must be between 0.0 < F < 1.0, with 1.0
- % corresponding to half the sample rate. They must be in increasing
- % order and start with 0.0 and end with 1.0.
- %
- % The filter B is real, and has linear phase, i.e., symmetric
- % coefficients obeying B(k) = B(N+2-k), k = 1,2,...,N+1.
- %
- % By default FIR2 windows the impulse response with a Hamming window.
- % Other available windows, including Boxcar, Hann, Bartlett, Blackman,
- % Kaiser and Chebwin can be specified with an optional trailing argument.
- % For example,
- % B = FIR2(N,F,A,bartlett(N+1)) uses a Bartlett window.
- % B = FIR2(N,F,A,chebwin(N+1,R)) uses a Chebyshev window.
- %
- % For filters with a gain other than zero at Fs/2, e.g., highpass
- % and bandstop filters, N must be even. Otherwise, N will be
- % incremented by one. In this case the window length should be
- % specified as N+2.
- %
- % See also FIR1, FIRLS, CFIRPM, FIRPM, BUTTER, CHEBY1, CHEBY2, YULEWALK,
- % FREQZ, FILTER.
- % Author(s): L. Shure, 3-27-87
- % C. Denham, 7-26-90, revised
- % Copyright 1988-2004 The MathWorks, Inc.
- % $Revision: 1.12.4.4 $ $Date: 2004/10/18 21:08:07 $
- % Optional input arguments (see the user's guide):
- % npt - number for interpolation
- % lap - the number of points for jumps in amplitude
- nargchk(3,6,nargin);
- [nn,msg1,msg2] = firchk(nn,ff(end),aa);
- error(msg1);
- warning(msg2);
- % Work with filter length instead of filter order
- nn = nn + 1;
- if (nargin > 3)
- if nargin == 4
- if length(npt) == 1
- if (2 ^ round(log(npt)/log(2))) ~= npt
- % NPT is not an even power of two
- npt = 2^ceil(log(npt)/log(2));
- end
- wind = hamming(nn);
- else
- wind = npt;
- if nn < 1024
- npt = 512;
- else
- npt = 2.^ceil(log(nn)/log(2));
- end
- end
- lap = fix(npt/25);
- elseif nargin == 5
- if length(npt) == 1
- if (2 ^ round(log(npt)/log(2))) ~= npt
- % NPT is not an even power of two
- npt = 2.^ceil(log(npt)/log(2));
- end
- if length(lap) == 1
- wind = hamming(nn);
- else
- wind = lap;
- lap = fix(npt/25);
- end
- else
- wind = npt;
- npt = lap;
- lap = fix(npt/25);
- end
- end
- elseif nargin == 3
- if nn < 1024
- npt = 512;
- else
- npt = 2.^ceil(log(nn)/log(2));
- end
- wind = hamming(nn);
- lap = fix(npt/25);
- end
- if nn ~= length(wind)
- error('The specified window must be the same as the filter length')
- end
- [mf,nf] = size(ff);
- [ma,na] = size(aa);
- if ma ~= mf || na ~= nf
- error('You must specify the same number of frequencies and amplitudes')
- end
- nbrk = max(mf,nf);
- if mf < nf
- ff = ff';
- aa = aa';
- end
- if abs(ff(1)) > eps || abs(ff(nbrk) - 1) > eps
- error('The first frequency must be 0 and the last 1')
- end
- % interpolate breakpoints onto large grid
- H = zeros(1,npt);
- nint=nbrk-1;
- df = diff(ff'); %#ok
- if (any(df < 0))
- error('Frequencies must be non-decreasing')
- end
- npt = npt + 1; % Length of [dc 1 2 ... nyquist] frequencies.
- nb = 1;
- H(1)=aa(1);
- for i=1:nint
- if df(i) == 0
- nb = nb - lap/2;
- ne = nb + lap;
- else
- ne = fix(ff(i+1)*npt);
- end
- if (nb < 0 || ne > npt)
- error('Too abrupt an amplitude change near end of frequency interval')
- end
- j=nb:ne;
- if nb == ne
- inc = 0;
- else
- inc = (j-nb)/(ne-nb);
- end
- H(nb:ne) = inc*aa(i+1) + (1 - inc)*aa(i);
- nb = ne + 1;
- end
- % Fourier time-shift.
- dt = 0.5 .* (nn - 1);
- rad = -dt .* sqrt(-1) .* pi .* (0:npt-1) ./ (npt-1);
- H = H .* exp(rad);
- H = [H conj(H(npt-1:-1:2))]; % Fourier transform of real series.
- ht = real(ifft(H)); % Symmetric real series.
- b = ht(1:nn); % Raw numerator.
- b = b .* wind(:).'; % Apply window.
- a = 1; % Denominator.
复制代码
[ 本帖最后由 sogooda 于 2008-4-29 21:52 编辑 ] |
评分
-
1
查看全部评分
-
|