1 ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
16 ## usage: b = fir1(n, w [, type] [, window] [, noscale])
18 ## Produce an order n FIR filter with the given frequency cutoff,
19 ## returning the n+1 filter coefficients in b.
21 ## n: order of the filter (1 less than the length of the filter)
23 ## strictly increasing vector in range [0, 1]
24 ## singleton for highpass or lowpass, vector pair for bandpass or
25 ## bandstop, or vector for alternating pass/stop filter.
26 ## type: choose between pass and stop bands
27 ## 'high' for highpass filter, cutoff at w
28 ## 'stop' for bandstop filter, edges at w = [lo, hi]
29 ## 'DC-0' for bandstop as first band of multiband filter
30 ## 'DC-1' for bandpass as first band of multiband filter
31 ## window: smoothing window
32 ## defaults to hamming(n+1) row vector
33 ## returned filter is the same shape as the smoothing window
34 ## noscale: choose whether to normalize or not
35 ## 'scale': set the magnitude of the center of the first passband to 1
36 ## 'noscale': don't normalize
38 ## To apply the filter, use the return vector b:
42 ## freqz(fir1(40,0.3));
43 ## freqz(fir1(15,[0.2, 0.5], 'stop')); # note the zero-crossing at 0.1
44 ## freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));
46 ## TODO: Consider using exact expression (in terms of sinc) for the
47 ## TODO: impulse response rather than relying on fir2.
48 ## TODO: Find reference to the requirement that order be even for
49 ## TODO: filters that end high. Figure out what to do with the
50 ## TODO: window in these cases
51 function b = fir1(n, w, varargin)
53 if nargin < 2 || nargin > 5
57 ## Assign default window, filter type and scale.
58 ## If single band edge, the first band defaults to a pass band to
59 ## create a lowpass filter. If multiple band edges, the first band
60 ## defaults to a stop band so that the two band case defaults to a
61 ## band pass filter. Ick.
64 ftype = (length(w)==1);
66 ## sort arglist, normalize any string
67 for i=1:length(varargin)
69 if ischar(arg), arg=lower(arg);end
70 if isempty(arg) continue; end # octave bug---can't switch on []
72 case {'low','stop','dc-1'}, ftype = 1;
73 case {'high','pass','bandpass','dc-0'}, ftype = 0;
74 case {'scale'}, scale = 1;
75 case {'noscale'}, scale = 0;
76 otherwise window = arg;
80 ## build response function according to fir2 requirements
83 f(1) = 0; f(2*bands)=1;
87 m(1:2:2*bands) = rem([1:bands]-(1-ftype),2);
88 m(2:2:2*bands) = m(1:2:2*bands);
90 ## Increment the order if the final band is a pass band. Something
91 ## about having a nyquist frequency of zero causing problems.
92 if rem(n,2)==1 && m(2*bands)==1,
93 warning("n must be even for highpass and bandstop filters. Incrementing.");
95 if isvector(window) && isreal(window) && !ischar(window)
96 ## Extend the window using interpolation
99 window = [window; window];
101 window = interp1(linspace(0,1,M),window,linspace(0,1,M+1),'linear');
103 window = interp1(linspace(0,1,M),window,linspace(0,1,M+1),'spline');
108 ## compute the filter
109 b = fir2(n, f, m, 512, 2, window);
111 ## normalize filter magnitude
113 ## find the middle of the first band edge
114 ## find the frequency of the normalizing gain
116 ## if the first band is a passband, use DC gain
119 ## for a highpass filter,
120 ## use the gain at half the sample frequency
123 ## otherwise, use the gain at the center
124 ## frequency of the first passband
125 w_o = f(3) + (f(4)-f(3))/2;
128 ## compute |h(w_o)|^-1
129 renorm = 1/abs(polyval(b, exp(-1i*pi*w_o)));
131 ## normalize the filter
137 %! freqz(fir1(40,0.3));
139 %! freqz(fir1(15,[0.2, 0.5], 'stop')); # note the zero-crossing at 0.1
141 %! freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));
143 %!assert(fir1(2, .5, 'low', @hanning, 'scale'), [0 1 0]');
144 %!assert(fir1(2, .5, 'low', "hanning", 'scale'), [0 1 0]');
145 %!assert(fir1(2, .5, 'low', hanning(3), 'scale'), [0 1 0]');
147 %!assert(fir1(10,.5,'noscale'), fir1(10,.5,'low','hamming','noscale'));
148 %!assert(fir1(10,.5,'high'), fir1(10,.5,'high','hamming','scale'));
149 %!assert(fir1(10,.5,'boxcar'), fir1(10,.5,'low','boxcar','scale'));
150 %!assert(fir1(10,.5,'hanning','scale'), fir1(10,.5,'scale','hanning','low'));
151 %!assert(fir1(10,.5,'haNNing','NOscale'), fir1(10,.5,'noscale','Hanning','LOW'));
152 %!assert(fir1(10,.5,'boxcar',[]), fir1(10,.5,'boxcar'));