|
0
|
1 ## Copyright (C) 2000 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2 of the License, or |
|
|
6 ## (at your option) any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, |
|
|
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 ## GNU General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; if not, write to the Free Software |
|
|
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
|
16 |
|
|
17 ## usage: b = fir1(n, w [, type] [, window] [, noscale]) |
|
|
18 ## |
|
|
19 ## Produce an order n FIR filter with the given frequency cutoff, |
|
|
20 ## returning the n+1 filter coefficients in b. |
|
|
21 ## |
|
|
22 ## n: order of the filter (1 less than the length of the filter) |
|
|
23 ## w: band edges |
|
|
24 ## strictly increasing vector in range [0, 1] |
|
|
25 ## singleton for highpass or lowpass, vector pair for bandpass or |
|
|
26 ## bandstop, or vector for alternating pass/stop filter. |
|
|
27 ## type: choose between pass and stop bands |
|
|
28 ## 'high' for highpass filter, cutoff at w |
|
|
29 ## 'stop' for bandstop filter, edges at w = [lo, hi] |
|
|
30 ## 'DC-0' for bandstop as first band of multiband filter |
|
|
31 ## 'DC-1' for bandpass as first band of multiband filter |
|
|
32 ## window: smoothing window |
|
|
33 ## defaults to hamming(n+1) row vector |
|
|
34 ## returned filter is the same shape as the smoothing window |
|
|
35 ## noscale: choose whether to normalize or not |
|
|
36 ## 'scale': set the magnitude of the center of the first passband to 1 |
|
|
37 ## 'noscale': don't normalize |
|
|
38 ## |
|
|
39 ## To apply the filter, use the return vector b: |
|
|
40 ## y=filter(b,1,x); |
|
|
41 ## |
|
|
42 ## Examples: |
|
|
43 ## freqz(fir1(40,0.3)); |
|
|
44 ## freqz(fir1(15,[0.2, 0.5], 'stop')); # note the zero-crossing at 0.1 |
|
|
45 ## freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale')); |
|
|
46 |
|
|
47 ## TODO: Consider using exact expression (in terms of sinc) for the |
|
|
48 ## TODO: impulse response rather than relying on fir2. |
|
|
49 ## TODO: Find reference to the requirement that order be even for |
|
|
50 ## TODO: filters that end high. Figure out what to do with the |
|
|
51 ## TODO: window in these cases---duplicating the central value |
|
|
52 ## TODO: like I currently do seems very reasonable. |
|
|
53 function b = fir1(n, w, ftype, window, scale) |
|
|
54 |
|
|
55 if nargin < 2 || nargin > 5 |
|
|
56 usage("b = fir1(n, w [, type] [, window] [, noscale])"); |
|
|
57 endif |
|
|
58 |
|
|
59 ## interpret arguments |
|
|
60 if nargin==2 |
|
|
61 ftype=[]; window=[]; scale=[]; |
|
|
62 elseif nargin==3 |
|
|
63 window=[]; scale=[]; |
|
|
64 if !isstr(ftype), window=ftype; ftype=[]; endif |
|
|
65 elseif nargin==4 |
|
|
66 scale=[]; |
|
|
67 if isstr(window), scale=window; window=[]; endif |
|
|
68 if !isstr(ftype), window=ftype; ftype=[]; endif |
|
|
69 endif |
|
|
70 |
|
|
71 ## If single band edge, the first band defaults to a pass band |
|
|
72 ## to create a lowpass filter. If multiple band edges, assume |
|
|
73 ## the first band is a stop band, so that the two band case defaults |
|
|
74 ## to a band pass filter. Ick. |
|
|
75 ftype = tolower(ftype); |
|
|
76 if isempty(ftype), ftype = length(w)==1; |
|
|
77 elseif strcmp(ftype, 'low'), ftype = 1; |
|
|
78 elseif strcmp(ftype, 'high'), ftype = 0; |
|
|
79 elseif strcmp(ftype, 'pass'), ftype = 0; |
|
|
80 elseif strcmp(ftype, 'stop'), ftype = 1; |
|
|
81 elseif strcmp(ftype, 'dc-0'), ftype = 0; |
|
|
82 elseif strcmp(ftype, 'dc-1'), ftype = 1; |
|
|
83 elseif isstr(ftype) |
|
|
84 error(["fir1 invalid filter type ", ftype]); |
|
|
85 else |
|
|
86 error("fir1 filter type should be a string"); |
|
|
87 endif |
|
|
88 |
|
|
89 ## scale the magnitude by default |
|
|
90 if isempty(scale) || strcmp(scale, 'scale'), scale = 1; |
|
|
91 elseif strcmp(scale, 'noscale'), scale=0; |
|
|
92 else error("fir1 scale must be 'scale' or 'noscale'"); |
|
|
93 endif |
|
|
94 |
|
|
95 ## use fir2 default filter |
|
|
96 if isempty(window) && isempty(scale), window = []; endif |
|
|
97 |
|
|
98 ## build response function according to fir2 requirements |
|
|
99 bands = length(w)+1; |
|
|
100 f = zeros(1,2*bands); |
|
|
101 f(1) = 0; f(2*bands)=1; |
|
|
102 f(2:2:2*bands-1) = w; |
|
|
103 f(3:2:2*bands-1) = w; |
|
|
104 m = zeros(1,2*bands); |
|
|
105 m(1:2:2*bands) = rem([1:bands]-(1-ftype),2); |
|
|
106 m(2:2:2*bands) = m(1:2:2*bands); |
|
|
107 |
|
|
108 ## Increment the order if the final band is a pass band. Something |
|
|
109 ## about having a nyquist frequency of zero causing problems. |
|
|
110 if rem(n,2)==1 && m(2*bands)==1, |
|
|
111 warning("n must be even for highpass and bandstop filters. Incrementing."); |
|
|
112 n=n+1; |
|
|
113 if !isempty(window), |
|
|
114 if rows(window) == 1 |
|
|
115 window = [window(1:n/2), window(n/2:n-1)]; |
|
|
116 else |
|
|
117 window = [window(1:n/2); window(n/2:n-1)]; |
|
|
118 endif |
|
|
119 endif |
|
|
120 endif |
|
|
121 |
|
|
122 ## compute the filter |
|
|
123 b = fir2(n, f, m, 512, 2, window); |
|
|
124 |
|
|
125 ## normalize filter magnitude |
|
|
126 if scale == 1 |
|
|
127 ## find the middle of the first band edge |
|
|
128 if m(1) == 1, w_o = (f(2)-f(1))/2; |
|
|
129 else w_o = f(3) + (f(4)-f(3))/2; |
|
|
130 endif |
|
|
131 |
|
|
132 ## compute |h(w_o)|^-1 |
|
|
133 renorm = 1/abs(polyval(b, exp(-1i*pi*w_o))); |
|
|
134 |
|
|
135 ## normalize the filter |
|
|
136 b = renorm*b; |
|
|
137 endif |
|
|
138 endfunction |
|
|
139 |
|
|
140 %!demo |
|
|
141 %! freqz(fir1(40,0.3)); |
|
|
142 %!demo |
|
|
143 %! freqz(fir1(15,[0.2, 0.5], 'stop')); # note the zero-crossing at 0.1 |
|
|
144 %!demo |
|
|
145 %! freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale')); |
