matlab程序中的fir1和fir2有什么不同

byrsky

最近看到一个matlab程序中同时出现fir1，fir2，但是我在help里面搜索有什么区别时，我却没发现二者有什么不同。谁能告诉我一下，fir1和fir2有什么不同，谢谢

[ 本帖最后由 eight 于 2008-4-29 21:10 编辑 ]

eight

原帖由 byrsky 于 2008-4-29 20:40 发表

                               
登录/注册后可看大图

最近看到一个matlab程序中同时出现fir1，fir2，但是我在help里面搜索有什么区别时，我却没发现二者有什么不同。谁能告诉我一下，fir1和fir2有什么不同，谢谢

这是自定义的变量或者函数吧？matlab 自带的库中应该没有 fir1、fir2 的。若确实如此，请附上代码

[ 本帖最后由 eight 于 2008-4-30 21:48 编辑 ]

sogooda

我的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

1. function [b,a] = fir1(N,Wn,varargin)
   
2. %FIR1   FIR filter design using the window method.
   
3. %   B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter
   
4. %   and returns the filter coefficients in length N+1 vector B.
   
5. %   The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0
   
6. %   corresponding to half the sample rate.  The filter B is real and
   
7. %   has linear phase.  The normalized gain of the filter at Wn is
   
8. %   -6 dB.
   
9. %
   
10. %   B = FIR1(N,Wn,'high') designs an N'th order highpass filter.
    
11. %   You can also use B = FIR1(N,Wn,'low') to design a lowpass filter.
    
12. %
    
13. %   If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an
    
14. %   order N bandpass filter with passband  W1 < W < W2. You can
    
15. %   also specify B = FIR1(N,Wn,'bandpass').  If Wn = [W1 W2],
    
16. %   B = FIR1(N,Wn,'stop') will design a bandstop filter.
    
17. %
    
18. %   If Wn is a multi-element vector,
    
19. %          Wn = [W1 W2 W3 W4 W5 ... WN],
    
20. %   FIR1 returns an order N multiband filter with bands
    
21. %    0 < W < W1, W1 < W < W2, ..., WN < W < 1.
    
22. %   B = FIR1(N,Wn,'DC-1') makes the first band a passband.
    
23. %   B = FIR1(N,Wn,'DC-0') makes the first band a stopband.
    
24. %
    
25. %   B = FIR1(N,Wn,WIN) designs an N-th order FIR filter using
    
26. %   the N+1 length vector WIN to window the impulse response.
    
27. %   If empty or omitted, FIR1 uses a Hamming window of length N+1.
    
28. %   For a complete list of available windows, see the help for the
    
29. %   WINDOW function. KAISER and CHEBWIN can be specified with an
    
30. %   optional trailing argument.  For example, B = FIR1(N,Wn,kaiser(N+1,4))
    
31. %   uses a Kaiser window with beta=4. B = FIR1(N,Wn,'high',chebwin(N+1,R))
    
32. %   uses a Chebyshev window with R decibels of relative sidelobe
    
33. %   attenuation.
    
34. %
    
35. %   For filters with a gain other than zero at Fs/2, e.g., highpass
    
36. %   and bandstop filters, N must be even.  Otherwise, N will be
    
37. %   incremented by one.  In this case the window length should be
    
38. %   specified as N+2.
    
39. %   
    
40. %   By default, the filter is scaled so the center of the first pass band
    
41. %   has magnitude exactly one after windowing. Use a trailing 'noscale'
    
42. %   argument to prevent this scaling, e.g. B = FIR1(N,Wn,'noscale'),
    
43. %   B = FIR1(N,Wn,'high','noscale'), B = FIR1(N,Wn,wind,'noscale').  You
    
44. %   can also specify the scaling explicitly, e.g. FIR1(N,Wn,'scale'), etc.
    
45. %
    
46. %   See also KAISERORD, FIRCLS1, FIR2, FIRLS, FIRCLS, CFIRPM,
    
47. %            FIRPM, FREQZ, FILTER, WINDOW.
    
48. 
    
49. %   FIR1 is an M-file implementation of program 5.2 in the IEEE
    
50. %   Programs for Digital Signal Processing tape.
    
51. 
    
52. %   Author(s): L. Shure
    
53. %              L. Shure, 4-5-90, revised
    
54. %              T. Krauss, 3-5-96, revised
    
55. %   Copyright 1988-2004 The MathWorks, Inc.
    
56. %   $Revision: 1.15.4.3 $  $Date: 2004/10/18 21:08:06 $
    
57. 
    
58. %   Reference(s):
    
59. %     [1] "Programs for Digital Signal Processing", IEEE Press
    
60. %         John Wiley & Sons, 1979, pg. 5.2-1.
    
61. 
    
62. error(nargchk(2,5,nargin));
    
63. 
    
64. % Parse optional input arguments
    
65. [Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin{:});
    
66. error(msg);
    
67. 
    
68. % Compute the frequency vector
    
69. [nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype);
    
70. error(msg);
    
71. 
    
72. % Compute the magnitude vector
    
73. [aa,First_Band] = desiredmag(Ftype,nbands);
    
74. 
    
75. % Check for appropriate filter order, increase when necessary
    
76. [N,msg1,msg2] = firchk(N,ff(end),aa);
    
77. error(msg1);
    
78. warning(msg2);
    
79. 
    
80. % Work with filter length (= order + 1)
    
81. L = N + 1;
    
82. 
    
83. % Check for valid window, or assign default if empty
    
84. [Wind,msg] = chkwindow(Wind,L);
    
85. error(msg);
    
86. 
    
87. % Compute unwindowed impulse response
    
88. hh = firls(L-1,ff,aa);
    
89. 
    
90. % Window impulse response to get the filter
    
91. b = hh.*Wind(:)';
    
92. a = 1;
    
93. 
    
94. if SCALING,   
    
95.     % Scale so that passband is approx 1
    
96.     b = scalefilter(b,First_Band,ff,L);   
    
97. end
    
98. 
    
99. %-------------------------------------------------------------------------
    
100. function [Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin)
     
101. %PARSEOPTARGS   Parse optional input arguments.
     
102. 
     
103. % Up to 3 optional input arguments, always in this order:
     
104. %   1 - Filter type flag, can be 'low','high','bandpass','stop','DC-0','DC-1'
     
105. %   2 - Window vector
     
106. %   3 - 'noscale' flag
     
107. 
     
108. % Initialize output args.
     
109. SCALING = [];
     
110. 
     
111. [Ftype,Wind,Scale] = assignoptargs(Wn,varargin{:});
     
112. 
     
113. [Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale);
     
114. if ~isempty(msg),
     
115.     return
     
116. end
     
117. 
     
118. switch lower(Scale),
     
119. case 'noscale';
     
120.     SCALING = 0;
     
121. case 'scale';
     
122.     SCALING = 1;
     
123. end
     
124. 
     
125. %--------------------------------------------------------------------------
     
126. function [Ftype,Wind,Scale] = assignoptargs(Wn,varargin)
     
127. %ASSIGNOPTARGS  Assign optional input arguments to the appropriate variables.
     
128. 
     
129. % default optional parameter values:
     
130. Wind = [];
     
131. Scale = 'scale';
     
132. Ftype = defaultftype(Wn);
     
133. 
     
134. 
     
135. switch length(varargin)
     
136. case 1
     
137.     if ischar(varargin{1}) && (length(varargin{1})>0),
     
138.         s = upper(varargin{1});
     
139.         switch upper(s)
     
140.         case {'SCALE','NOSCALE'}
     
141.             Scale = s;
     
142.         otherwise
     
143.             Ftype = s;
     
144.         end
     
145.     else
     
146.         Wind = varargin{1};
     
147.     end
     
148. case 2
     
149.     if ischar(varargin{1})
     
150.         Ftype = varargin{1};
     
151.     else
     
152.         Wind = varargin{1};
     
153.     end
     
154.     if ischar(varargin{2})
     
155.         Scale = varargin{2};
     
156.     else
     
157.         Wind = varargin{2};
     
158.     end
     
159. case 3
     
160.     Ftype = varargin{1};
     
161.     Wind = varargin{2};
     
162.     Scale = varargin{3};
     
163. end
     
164. 
     
165. %--------------------------------------------------------------------------
     
166. function [Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale)
     
167. %VALIDATEARGS  Test if arguments are valid.
     
168. 
     
169. msg = '';
     
170. 
     
171. % Assign a default Ftype when an empty is given. Backwards compatibility
     
172. if isempty(Ftype),
     
173.     Ftype = defaultftype(Wn);
     
174. end
     
175. 
     
176. 
     
177. Ftypeopts = {'LOW','HIGH','BANDPASS','STOP','DC-0','DC-1'};
     
178. Scaleopts = {'NOSCALE','SCALE'};
     
179. 
     
180. 
     
181. indx = find(strncmpi(Ftype, Ftypeopts, length(Ftype)));
     
182. if isempty(indx),
     
183.     msg = 'Unrecognized or ambiguous filter type specified.';
     
184.     return
     
185. else
     
186.     Ftype = Ftypeopts{indx};
     
187. end
     
188. 
     
189. scaleindx = find(strncmpi(Scale, Scaleopts, length(Scale)));
     
190. if isempty(scaleindx),
     
191.     msg = 'Scaling option must be ''noscale'' or ''scale''.';
     
192.     return
     
193. else
     
194.     Scale = Scaleopts{scaleindx};
     
195. end
     
196. 
     
197. if ~any(size(Wind) <= 1),
     
198.     msg = 'The window specified must be a vector';
     
199.     return
     
200. else
     
201.     Wind = Wind(:).'; % Make it a row vector
     
202. end
     
203. 
     
204. %--------------------------------------------------------------------------
     
205. function [nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype)
     
206. %DESIREDFREQ  Compute the vector of frequencies to pass to FIRLS.
     
207. %
     
208. %   Inputs:
     
209. %           Wn    - vector of cutoff frequencies.
     
210. %           Ftype - string with desired response ('low','high',...)
     
211. %
     
212. %   Outputs:
     
213. %           nbands - number of frequency bands.
     
214. %           ff     - vector of frequencies to pass to FIRLS.
     
215. %           Ftype  - converted filter type (if it's necessary to convert)
     
216. 
     
217. % Initialize output args.
     
218. nbands = [];
     
219. ff     = [];
     
220. msg    = '';
     
221. 
     
222. 
     
223. if  any( Wn<0 | Wn>1 ),
     
224.    msg = 'Frequencies must fall in range between 0 and 1.';
     
225.    return
     
226. end
     
227. if  any(diff(Wn)<0),
     
228.    msg = 'Frequencies must be increasing';
     
229.    return
     
230. end
     
231. 
     
232. Wn = Wn(:)';
     
233. 
     
234. nbands = length(Wn) + 1;
     
235. 
     
236. if (nbands > 2) && strcmpi(Ftype,'bandpass'),
     
237.     Ftype = 'DC-0';  % make sure default 3 band filter is bandpass
     
238. end
     
239. 
     
240. ff = [0,Wn(1:nbands-1); Wn(1:nbands-1),1];
     
241. 
     
242. ff = ff(:);
     
243. 
     
244. %-------------------------------------------------------------------------
     
245. function [aa,First_Band] = desiredmag(Ftype,nbands)
     
246. %DESIREDMAG  Compute the magnitude vector to pass to FIRLS.
     
247. 
     
248. First_Band = isempty(findstr('DC-0',Ftype)) && isempty(findstr('HIGH',Ftype));
     
249. mags = rem( First_Band + (0:nbands-1), 2);
     
250. aa = [mags(:)'; mags(:)'];
     
251. 
     
252. aa = aa(:);
     
253. %--------------------------------------------------------------------------
     
254. function [Wind,msg] = chkwindow(Wind,L)
     
255. %CHKWINDOW   Check if specified window is valid, assign default if empty.
     
256. 
     
257. msg = '';
     
258. 
     
259. if isempty(Wind),
     
260.    % Replace the following with the default window of your choice.
     
261.    Wind = hamming(L);
     
262. end
     
263.    
     
264. if length(Wind) ~= L
     
265.    msg = 'The window length must be the same as the filter length.';
     
266. end
     
267. %
     
268. % to use Kaiser window, beta must be supplied
     
269. % att = 60; % dB of attenuation desired in sidelobe
     
270. % beta = 0.1102*(att-8.7);
     
271. % wind = kaiser(L,beta);
     
272. 
     
273. %---------------------------------------------------------------------------
     
274. function b = scalefilter(b,First_Band,ff,L)
     
275. %SCALEFILTER   Scale fitler to have passband approx. equal to one.
     
276. 
     
277. if First_Band
     
278.     b = b / sum(b);  % unity gain at DC
     
279. else
     
280.     if ff(4)==1
     
281.         % unity gain at Fs/2
     
282.         f0 = 1;
     
283.     else
     
284.         % unity gain at center of first passband
     
285.         f0 = mean(ff(3:4));
     
286.     end
     
287.     b = b / abs( exp(-j*2*pi*(0:L-1)*(f0/2))*(b.') );
     
288. end
     
289. 
     
290. %----------------------------------------------------------------------------
     
291. function Ftype = defaultftype(Wn)
     
292. %DEFAULTFTYPE  Assign default filter type depending on number of bands.
     
293. 
     
294. if length(Wn) == 1,
     
295.     Ftype = 'low';
     
296. elseif length(Wn) == 2,
     
297.     Ftype = 'bandpass';
     
298. elseif length(Wn) >= 3,
     
299.     Ftype = 'dc-0';
     
300. end

复制代码

fir2.m

1. function [b,a] = fir2(nn, ff, aa, npt, lap, wind)
   
2. %FIR2   FIR arbitrary shape filter design using the frequency sampling method.
   
3. %   B = FIR2(N,F,A) designs an N'th order FIR digital filter with the
   
4. %   frequency response specified by vectors F and A, and returns the
   
5. %   filter coefficients in length N+1 vector B.  Vectors F and A specify
   
6. %   the frequency and magnitude breakpoints for the filter such that
   
7. %   PLOT(F,A) would show a plot of the desired frequency response.
   
8. %   The frequencies in F must be between 0.0 < F < 1.0, with 1.0
   
9. %   corresponding to half the sample rate. They must be in increasing
   
10. %   order and start with 0.0 and end with 1.0.
    
11. %
    
12. %   The filter B is real, and has linear phase, i.e., symmetric
    
13. %   coefficients obeying B(k) =  B(N+2-k), k = 1,2,...,N+1.
    
14. %
    
15. %   By default FIR2 windows the impulse response with a Hamming window.
    
16. %   Other available windows, including Boxcar, Hann, Bartlett, Blackman,
    
17. %   Kaiser and Chebwin can be specified with an optional trailing argument.
    
18. %   For example,
    
19. %   B = FIR2(N,F,A,bartlett(N+1)) uses a Bartlett window.
    
20. %   B = FIR2(N,F,A,chebwin(N+1,R)) uses a Chebyshev window.
    
21. %
    
22. %   For filters with a gain other than zero at Fs/2, e.g., highpass
    
23. %   and bandstop filters, N must be even.  Otherwise, N will be
    
24. %   incremented by one.  In this case the window length should be
    
25. %   specified as N+2.
    
26. %
    
27. %   See also FIR1, FIRLS, CFIRPM, FIRPM, BUTTER, CHEBY1, CHEBY2, YULEWALK,
    
28. %   FREQZ, FILTER.
    
29. 
    
30. %   Author(s): L. Shure, 3-27-87
    
31. %         C. Denham, 7-26-90, revised
    
32. %   Copyright 1988-2004 The MathWorks, Inc.
    
33. %   $Revision: 1.12.4.4 $  $Date: 2004/10/18 21:08:07 $
    
34. 
    
35. %   Optional input arguments (see the user's guide):
    
36. %     npt - number for interpolation
    
37. %     lap - the number of points for jumps in amplitude
    
38. 
    
39. 
    
40. nargchk(3,6,nargin);
    
41. 
    
42. [nn,msg1,msg2] = firchk(nn,ff(end),aa);
    
43. error(msg1);
    
44. warning(msg2);
    
45. 
    
46. % Work with filter length instead of filter order
    
47. nn = nn + 1;
    
48. 
    
49. if (nargin > 3)
    
50.     if nargin == 4
    
51.         if length(npt) == 1
    
52.             if (2 ^ round(log(npt)/log(2))) ~= npt
    
53.                 % NPT is not an even power of two
    
54.                 npt = 2^ceil(log(npt)/log(2));
    
55.             end
    
56.             wind = hamming(nn);
    
57.         else
    
58.             wind = npt;
    
59.             if nn < 1024
    
60.                 npt = 512;
    
61.             else
    
62.                 npt = 2.^ceil(log(nn)/log(2));
    
63.             end
    
64.         end
    
65.         lap = fix(npt/25);
    
66.     elseif nargin == 5
    
67.         if length(npt) == 1
    
68.             if (2 ^ round(log(npt)/log(2))) ~= npt
    
69.                 % NPT is not an even power of two
    
70.                 npt = 2.^ceil(log(npt)/log(2));
    
71.             end
    
72.             if length(lap) == 1
    
73.                 wind = hamming(nn);
    
74.             else
    
75.                 wind = lap;
    
76.                 lap = fix(npt/25);
    
77.             end
    
78.         else
    
79.             wind = npt;
    
80.             npt = lap;
    
81.             lap = fix(npt/25);
    
82.         end
    
83.     end
    
84. elseif nargin == 3
    
85.     if nn < 1024
    
86.         npt = 512;
    
87.     else
    
88.         npt = 2.^ceil(log(nn)/log(2));
    
89.     end
    
90.     wind = hamming(nn);
    
91.     lap = fix(npt/25);
    
92. end
    
93. 
    
94. if nn ~= length(wind)
    
95.     error('The specified window must be the same as the filter length')
    
96. end
    
97. 
    
98. [mf,nf] = size(ff);
    
99. [ma,na] = size(aa);
    
100. if ma ~= mf || na ~= nf
     
101.     error('You must specify the same number of frequencies and amplitudes')
     
102. end
     
103. nbrk = max(mf,nf);
     
104. if mf < nf
     
105.     ff = ff';
     
106.     aa = aa';
     
107. end
     
108. 
     
109. if abs(ff(1)) > eps || abs(ff(nbrk) - 1) > eps
     
110.     error('The first frequency must be 0 and the last 1')
     
111. end
     
112. 
     
113. % interpolate breakpoints onto large grid
     
114. 
     
115. H = zeros(1,npt);
     
116. nint=nbrk-1;
     
117. df = diff(ff'); %#ok
     
118. if (any(df < 0))
     
119.     error('Frequencies must be non-decreasing')
     
120. end
     
121. 
     
122. npt = npt + 1;   % Length of [dc 1 2 ... nyquist] frequencies.
     
123. 
     
124. nb = 1;
     
125. H(1)=aa(1);
     
126. for i=1:nint
     
127.     if df(i) == 0
     
128.         nb = nb - lap/2;
     
129.         ne = nb + lap;
     
130.     else
     
131.         ne = fix(ff(i+1)*npt);
     
132.     end
     
133.     if (nb < 0 || ne > npt)
     
134.         error('Too abrupt an amplitude change near end of frequency interval')
     
135.     end
     
136.     j=nb:ne;
     
137.     if nb == ne
     
138.         inc = 0;
     
139.     else
     
140.         inc = (j-nb)/(ne-nb);
     
141.     end
     
142.     H(nb:ne) = inc*aa(i+1) + (1 - inc)*aa(i);
     
143.     nb = ne + 1;
     
144. end
     
145. 
     
146. % Fourier time-shift.
     
147. 
     
148. dt = 0.5 .* (nn - 1);
     
149. rad = -dt .* sqrt(-1) .* pi .* (0:npt-1) ./ (npt-1);
     
150. H = H .* exp(rad);
     
151. 
     
152. H = [H conj(H(npt-1:-1:2))];   % Fourier transform of real series.
     
153. ht = real(ifft(H));            % Symmetric real series.
     
154. 
     
155. b = ht(1:nn);         % Raw numerator.
     
156. b = b .* wind(:).';   % Apply window.
     
157. 
     
158. a = 1;                % Denominator.

复制代码

[ 本帖最后由 sogooda 于 2008-4-29 21:52 编辑 ]

eight

原帖由 sogooda 于 2008-4-29 21:45 发表

                               
登录/注册后可看大图

我的matlabR14（V7.04）里确实有fir1和fir2，在Signal Processing Toolbox里。eight院长没有可能是matlab的版本问题吧。另外，楼主，关于二者的区别就在help文档的Note里写的很清楚了啊。

Note    Use fir2 for w ...

呵呵，最近比较忙，所以没有时间打开matlab查证:lol

byrsky

前辈可不可以用中文简单的介绍一下这两个函数的区别呀，英文看不大懂:@L

byrsky

谢谢院长同志帮我解决问题:@)

sigma665

看不懂可以查字典吧:@L