home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / FIR1.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.8 KB  |  156 lines
  1. function [b,a] = fir1(N,Wn,Ftype,Wind)
  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.
  7. %
  8. %    If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an
  9. %    order N bandpass filter with passband  W1 < W < W2.
  10. %    B = FIR1(N,Wn,'high') designs a highpass filter.
  11. %    B = FIR1(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
  12. %    For highpass and bandstop filters, N must be even.
  13. %    
  14. %    By default FIR1 uses a Hamming window.  Other available windows,
  15. %    including Boxcar, Hanning, Bartlett, Blackman, Kaiser and Chebwin
  16. %    can be specified with an optional trailing argument.  For example,
  17. %    B = FIR1(N,Wn,bartlett(N+1)) uses a Bartlett window.
  18. %    B = FIR1(N,Wn,'high',chebwin(N+1,R)) uses a Chebyshev window.
  19. %
  20. %    FIR1 is an M-file implementation of program 5.2 in the IEEE
  21. %    Programs for Digital Signal Processing tape.  See also FIR2,
  22. %    BUTTER, CHEBY1, CHEBY2, YULEWALK, FREQZ and FILTER.
  23.  
  24. %    References:
  25. %      "Programs for Digital Signal Processing", IEEE Press
  26. %      John Wiley & Sons, 1979, pg. 5.2-1.
  27.  
  28. %    Author: L. Shure
  29. %    Revised 4-MAY-90, LS.
  30. %    (c) Copyright 1987-90, by The MathWorks, Inc.
  31.  
  32. nw = 0;
  33. a = 1;
  34. if nargin == 3
  35.     if ~isstr(Ftype);
  36.         nw = max(size(Ftype));
  37.               Wind = Ftype;
  38.         Ftype = [];
  39.     end
  40. elseif nargin == 4
  41.    nw = max(size(Wind));
  42. else
  43.    Ftype = [];
  44. end
  45.  
  46. Btype = 1;
  47. if nargin > 2 & max(size(Ftype)) > 0
  48.     Btype = 3;
  49. end
  50. if max(size(Wn)) == 2
  51.     Btype = Btype + 1;
  52. end
  53.  
  54. N = N + 1;
  55. odd = rem(N, 2);
  56. if (Btype == 3 | Btype == 4)
  57.    if (~odd)
  58.       disp('For highpass and bandstop filters, order must be even.')
  59.       disp('Order is being increased by 1.')
  60.       N = N + 1;
  61.       odd = 1;
  62.    end
  63. end
  64. if nw ~= 0 & nw ~= N
  65.    error('The window length must be the same as the filter length.')
  66. end
  67. if nw > 0
  68.    wind = Wind;
  69. else             % replace the following with the default window of your choice.
  70.    wind = hamming(N);
  71. end
  72. %
  73. % to use Chebyshev window, you must specify ripple
  74. % ripple=60;
  75. % wind=chebwin(N, ripple);
  76. %
  77. % to use Kaiser window, beta must be supplied
  78. % att = 60; % dB of attenuation desired in sidelobe
  79. % beta = .1102*(att-8.7);
  80. % wind = kaiser(N,beta);
  81.  
  82. fl = Wn(1)/2;
  83. if (Btype == 2 | Btype == 4)
  84.    fh = Wn(2)/2;
  85.    if (fl >= .5 | fl <= 0 | fh >= .5 | fh <= 0.)
  86.       error('Frequencies must fall in range between 0 and 1.')
  87.    end
  88.    c1=fh-fl;
  89.    if (c1 <= 0)
  90.       error('Wn(1) must be less than Wn(2).')
  91.    end
  92. else
  93.    c1=fl;
  94.    if (fl >= .5 | fl <= 0)
  95.       error('Frequency must lie between 0 and 1')
  96.    end 
  97. end
  98.  
  99. nhlf = fix((N + 1)/2);
  100. i1=1 + odd;
  101.  
  102. if Btype == 1            % lowpass
  103. if odd
  104.    b(1) = 2*c1;
  105. end
  106. xn=(odd:nhlf-1) + .5*(1-odd);
  107. c=pi*xn;
  108. c3=2*c1*c;
  109. b(i1:nhlf)=(sin(c3)./c);
  110. b = real([b(nhlf:-1:i1) b(1:nhlf)].*wind');
  111. gain = abs(polyval(b,1));
  112. b = b/gain;
  113. return;
  114.  
  115. elseif Btype ==2        % bandpass
  116. b(1) = 2*c1;
  117. xn=(odd:nhlf-1)+.5*(1-odd);
  118. c=pi*xn;
  119. c3=c*c1;
  120. b(i1:nhlf)=2*sin(c3).*cos(c*(fl+fh))./c;
  121. b=real([b(nhlf:-1:i1) b(1:nhlf)].*wind');
  122. gain = abs(polyval(b,exp(sqrt(-1)*pi*(fl+fh))));
  123. b = b/gain;
  124. return;
  125.  
  126. elseif Btype == 3        % highpass
  127. b(1)=2*c1;
  128. xn=(odd:nhlf-1);
  129. c=pi*xn;
  130. c3=2*c1*c;
  131. b(i1:nhlf)=sin(c3)./c;
  132. att=60.;
  133. b=real([b(nhlf:-1:i1) b(1:nhlf)].*wind');
  134. b(nhlf)=1-b(nhlf);
  135. b(1:nhlf-1)=-b(1:nhlf-1);
  136. b(nhlf+1:N)=-b(nhlf+1:N);
  137. gain = abs(polyval(b,-1));
  138. b = b/gain;
  139. return;
  140.  
  141. elseif Btype == 4        % bandstop
  142. b(1) = 2*c1;
  143. xn=(odd:nhlf-1)+.5*(1-odd);
  144. c=pi*xn;
  145. c3=c*c1;
  146. b(i1:nhlf)=2*sin(c3).*cos(c*(fl+fh))./c;
  147. b=real([b(nhlf:-1:i1) b(1:nhlf)].*wind');
  148. b(nhlf)=1-b(nhlf);
  149. b(1:nhlf-1)=-b(1:nhlf-1);
  150. b(nhlf+1:N)=-b(nhlf+1:N);
  151. gain = abs(polyval(b,1));
  152. b = b/gain;
  153. return;
  154. end
  155.  
  156.