声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 11675|回复: 6

[综合讨论] matlab程序中的fir1和fir2有什么不同

[复制链接]
发表于 2008-4-29 20:40 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

[ 本帖最后由 eight 于 2008-4-29 21:10 编辑 ]
回复
分享到:

使用道具 举报

发表于 2008-4-29 21:11 | 显示全部楼层
原帖由 byrsky 于 2008-4-29 20:40 发表
最近看到一个matlab程序中同时出现fir1,fir2,但是我在help里面搜索有什么区别时,我却没发现二者有什么不同。谁能告诉我一下,fir1和fir2有什么不同,谢谢
这是自定义的变量或者函数吧?matlab 自带的库中应该没有 fir1、fir2 的。若确实如此,请附上代码

[ 本帖最后由 eight 于 2008-4-30 21:48 编辑 ]
发表于 2008-4-29 21:45 | 显示全部楼层

回复 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
  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. %   FIR1 is an M-file implementation of program 5.2 in the IEEE
  49. %   Programs for Digital Signal Processing tape.

  50. %   Author(s): L. Shure
  51. %              L. Shure, 4-5-90, revised
  52. %              T. Krauss, 3-5-96, revised
  53. %   Copyright 1988-2004 The MathWorks, Inc.
  54. %   $Revision: 1.15.4.3 $  $Date: 2004/10/18 21:08:06 $

  55. %   Reference(s):
  56. %     [1] "Programs for Digital Signal Processing", IEEE Press
  57. %         John Wiley & Sons, 1979, pg. 5.2-1.

  58. error(nargchk(2,5,nargin));

  59. % Parse optional input arguments
  60. [Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin{:});
  61. error(msg);

  62. % Compute the frequency vector
  63. [nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype);
  64. error(msg);

  65. % Compute the magnitude vector
  66. [aa,First_Band] = desiredmag(Ftype,nbands);

  67. % Check for appropriate filter order, increase when necessary
  68. [N,msg1,msg2] = firchk(N,ff(end),aa);
  69. error(msg1);
  70. warning(msg2);

  71. % Work with filter length (= order + 1)
  72. L = N + 1;

  73. % Check for valid window, or assign default if empty
  74. [Wind,msg] = chkwindow(Wind,L);
  75. error(msg);

  76. % Compute unwindowed impulse response
  77. hh = firls(L-1,ff,aa);

  78. % Window impulse response to get the filter
  79. b = hh.*Wind(:)';
  80. a = 1;

  81. if SCALING,   
  82.     % Scale so that passband is approx 1
  83.     b = scalefilter(b,First_Band,ff,L);   
  84. end

  85. %-------------------------------------------------------------------------
  86. function [Ftype,Wind,SCALING,msg] = parseoptargs(Wn,varargin)
  87. %PARSEOPTARGS   Parse optional input arguments.

  88. % Up to 3 optional input arguments, always in this order:
  89. %   1 - Filter type flag, can be 'low','high','bandpass','stop','DC-0','DC-1'
  90. %   2 - Window vector
  91. %   3 - 'noscale' flag

  92. % Initialize output args.
  93. SCALING = [];

  94. [Ftype,Wind,Scale] = assignoptargs(Wn,varargin{:});

  95. [Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale);
  96. if ~isempty(msg),
  97.     return
  98. end

  99. switch lower(Scale),
  100. case 'noscale';
  101.     SCALING = 0;
  102. case 'scale';
  103.     SCALING = 1;
  104. end

  105. %--------------------------------------------------------------------------
  106. function [Ftype,Wind,Scale] = assignoptargs(Wn,varargin)
  107. %ASSIGNOPTARGS  Assign optional input arguments to the appropriate variables.

  108. % default optional parameter values:
  109. Wind = [];
  110. Scale = 'scale';
  111. Ftype = defaultftype(Wn);


  112. switch length(varargin)
  113. case 1
  114.     if ischar(varargin{1}) && (length(varargin{1})>0),
  115.         s = upper(varargin{1});
  116.         switch upper(s)
  117.         case {'SCALE','NOSCALE'}
  118.             Scale = s;
  119.         otherwise
  120.             Ftype = s;
  121.         end
  122.     else
  123.         Wind = varargin{1};
  124.     end
  125. case 2
  126.     if ischar(varargin{1})
  127.         Ftype = varargin{1};
  128.     else
  129.         Wind = varargin{1};
  130.     end
  131.     if ischar(varargin{2})
  132.         Scale = varargin{2};
  133.     else
  134.         Wind = varargin{2};
  135.     end
  136. case 3
  137.     Ftype = varargin{1};
  138.     Wind = varargin{2};
  139.     Scale = varargin{3};
  140. end

  141. %--------------------------------------------------------------------------
  142. function [Ftype,Wind,Scale,msg] = validateargs(Wn,Ftype,Wind,Scale)
  143. %VALIDATEARGS  Test if arguments are valid.

  144. msg = '';

  145. % Assign a default Ftype when an empty is given. Backwards compatibility
  146. if isempty(Ftype),
  147.     Ftype = defaultftype(Wn);
  148. end


  149. Ftypeopts = {'LOW','HIGH','BANDPASS','STOP','DC-0','DC-1'};
  150. Scaleopts = {'NOSCALE','SCALE'};


  151. indx = find(strncmpi(Ftype, Ftypeopts, length(Ftype)));
  152. if isempty(indx),
  153.     msg = 'Unrecognized or ambiguous filter type specified.';
  154.     return
  155. else
  156.     Ftype = Ftypeopts{indx};
  157. end

  158. scaleindx = find(strncmpi(Scale, Scaleopts, length(Scale)));
  159. if isempty(scaleindx),
  160.     msg = 'Scaling option must be ''noscale'' or ''scale''.';
  161.     return
  162. else
  163.     Scale = Scaleopts{scaleindx};
  164. end

  165. if ~any(size(Wind) <= 1),
  166.     msg = 'The window specified must be a vector';
  167.     return
  168. else
  169.     Wind = Wind(:).'; % Make it a row vector
  170. end

  171. %--------------------------------------------------------------------------
  172. function [nbands,ff,Ftype,msg] = desiredfreq(Wn,Ftype)
  173. %DESIREDFREQ  Compute the vector of frequencies to pass to FIRLS.
  174. %
  175. %   Inputs:
  176. %           Wn    - vector of cutoff frequencies.
  177. %           Ftype - string with desired response ('low','high',...)
  178. %
  179. %   Outputs:
  180. %           nbands - number of frequency bands.
  181. %           ff     - vector of frequencies to pass to FIRLS.
  182. %           Ftype  - converted filter type (if it's necessary to convert)

  183. % Initialize output args.
  184. nbands = [];
  185. ff     = [];
  186. msg    = '';


  187. if  any( Wn<0 | Wn>1 ),
  188.    msg = 'Frequencies must fall in range between 0 and 1.';
  189.    return
  190. end
  191. if  any(diff(Wn)<0),
  192.    msg = 'Frequencies must be increasing';
  193.    return
  194. end

  195. Wn = Wn(:)';

  196. nbands = length(Wn) + 1;

  197. if (nbands > 2) && strcmpi(Ftype,'bandpass'),
  198.     Ftype = 'DC-0';  % make sure default 3 band filter is bandpass
  199. end

  200. ff = [0,Wn(1:nbands-1); Wn(1:nbands-1),1];

  201. ff = ff(:);

  202. %-------------------------------------------------------------------------
  203. function [aa,First_Band] = desiredmag(Ftype,nbands)
  204. %DESIREDMAG  Compute the magnitude vector to pass to FIRLS.

  205. First_Band = isempty(findstr('DC-0',Ftype)) && isempty(findstr('HIGH',Ftype));
  206. mags = rem( First_Band + (0:nbands-1), 2);
  207. aa = [mags(:)'; mags(:)'];

  208. aa = aa(:);
  209. %--------------------------------------------------------------------------
  210. function [Wind,msg] = chkwindow(Wind,L)
  211. %CHKWINDOW   Check if specified window is valid, assign default if empty.

  212. msg = '';

  213. if isempty(Wind),
  214.    % Replace the following with the default window of your choice.
  215.    Wind = hamming(L);
  216. end
  217.    
  218. if length(Wind) ~= L
  219.    msg = 'The window length must be the same as the filter length.';
  220. end
  221. %
  222. % to use Kaiser window, beta must be supplied
  223. % att = 60; % dB of attenuation desired in sidelobe
  224. % beta = 0.1102*(att-8.7);
  225. % wind = kaiser(L,beta);

  226. %---------------------------------------------------------------------------
  227. function b = scalefilter(b,First_Band,ff,L)
  228. %SCALEFILTER   Scale fitler to have passband approx. equal to one.

  229. if First_Band
  230.     b = b / sum(b);  % unity gain at DC
  231. else
  232.     if ff(4)==1
  233.         % unity gain at Fs/2
  234.         f0 = 1;
  235.     else
  236.         % unity gain at center of first passband
  237.         f0 = mean(ff(3:4));
  238.     end
  239.     b = b / abs( exp(-j*2*pi*(0:L-1)*(f0/2))*(b.') );
  240. end

  241. %----------------------------------------------------------------------------
  242. function Ftype = defaultftype(Wn)
  243. %DEFAULTFTYPE  Assign default filter type depending on number of bands.

  244. if length(Wn) == 1,
  245.     Ftype = 'low';
  246. elseif length(Wn) == 2,
  247.     Ftype = 'bandpass';
  248. elseif length(Wn) >= 3,
  249.     Ftype = 'dc-0';
  250. 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. %   Author(s): L. Shure, 3-27-87
  30. %         C. Denham, 7-26-90, revised
  31. %   Copyright 1988-2004 The MathWorks, Inc.
  32. %   $Revision: 1.12.4.4 $  $Date: 2004/10/18 21:08:07 $

  33. %   Optional input arguments (see the user's guide):
  34. %     npt - number for interpolation
  35. %     lap - the number of points for jumps in amplitude


  36. nargchk(3,6,nargin);

  37. [nn,msg1,msg2] = firchk(nn,ff(end),aa);
  38. error(msg1);
  39. warning(msg2);

  40. % Work with filter length instead of filter order
  41. nn = nn + 1;

  42. if (nargin > 3)
  43.     if nargin == 4
  44.         if length(npt) == 1
  45.             if (2 ^ round(log(npt)/log(2))) ~= npt
  46.                 % NPT is not an even power of two
  47.                 npt = 2^ceil(log(npt)/log(2));
  48.             end
  49.             wind = hamming(nn);
  50.         else
  51.             wind = npt;
  52.             if nn < 1024
  53.                 npt = 512;
  54.             else
  55.                 npt = 2.^ceil(log(nn)/log(2));
  56.             end
  57.         end
  58.         lap = fix(npt/25);
  59.     elseif nargin == 5
  60.         if length(npt) == 1
  61.             if (2 ^ round(log(npt)/log(2))) ~= npt
  62.                 % NPT is not an even power of two
  63.                 npt = 2.^ceil(log(npt)/log(2));
  64.             end
  65.             if length(lap) == 1
  66.                 wind = hamming(nn);
  67.             else
  68.                 wind = lap;
  69.                 lap = fix(npt/25);
  70.             end
  71.         else
  72.             wind = npt;
  73.             npt = lap;
  74.             lap = fix(npt/25);
  75.         end
  76.     end
  77. elseif nargin == 3
  78.     if nn < 1024
  79.         npt = 512;
  80.     else
  81.         npt = 2.^ceil(log(nn)/log(2));
  82.     end
  83.     wind = hamming(nn);
  84.     lap = fix(npt/25);
  85. end

  86. if nn ~= length(wind)
  87.     error('The specified window must be the same as the filter length')
  88. end

  89. [mf,nf] = size(ff);
  90. [ma,na] = size(aa);
  91. if ma ~= mf || na ~= nf
  92.     error('You must specify the same number of frequencies and amplitudes')
  93. end
  94. nbrk = max(mf,nf);
  95. if mf < nf
  96.     ff = ff';
  97.     aa = aa';
  98. end

  99. if abs(ff(1)) > eps || abs(ff(nbrk) - 1) > eps
  100.     error('The first frequency must be 0 and the last 1')
  101. end

  102. % interpolate breakpoints onto large grid

  103. H = zeros(1,npt);
  104. nint=nbrk-1;
  105. df = diff(ff'); %#ok
  106. if (any(df < 0))
  107.     error('Frequencies must be non-decreasing')
  108. end

  109. npt = npt + 1;   % Length of [dc 1 2 ... nyquist] frequencies.

  110. nb = 1;
  111. H(1)=aa(1);
  112. for i=1:nint
  113.     if df(i) == 0
  114.         nb = nb - lap/2;
  115.         ne = nb + lap;
  116.     else
  117.         ne = fix(ff(i+1)*npt);
  118.     end
  119.     if (nb < 0 || ne > npt)
  120.         error('Too abrupt an amplitude change near end of frequency interval')
  121.     end
  122.     j=nb:ne;
  123.     if nb == ne
  124.         inc = 0;
  125.     else
  126.         inc = (j-nb)/(ne-nb);
  127.     end
  128.     H(nb:ne) = inc*aa(i+1) + (1 - inc)*aa(i);
  129.     nb = ne + 1;
  130. end

  131. % Fourier time-shift.

  132. dt = 0.5 .* (nn - 1);
  133. rad = -dt .* sqrt(-1) .* pi .* (0:npt-1) ./ (npt-1);
  134. H = H .* exp(rad);

  135. H = [H conj(H(npt-1:-1:2))];   % Fourier transform of real series.
  136. ht = real(ifft(H));            % Symmetric real series.

  137. b = ht(1:nn);         % Raw numerator.
  138. b = b .* wind(:).';   % Apply window.

  139. a = 1;                % Denominator.
复制代码

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

评分

1

查看全部评分

发表于 2008-4-30 11:04 | 显示全部楼层
原帖由 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
 楼主| 发表于 2008-4-30 20:12 | 显示全部楼层

回复 3楼 的帖子

前辈可不可以用中文简单的介绍一下这两个函数的区别呀,英文看不大懂:@L
 楼主| 发表于 2008-4-30 20:13 | 显示全部楼层

回复 2楼 的帖子

谢谢院长同志帮我解决问题:@)
发表于 2008-4-30 21:48 | 显示全部楼层

回复 5楼 的帖子

看不懂可以查字典吧:@L
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 02:48 , Processed in 0.095305 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表