home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / FIR2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.6 KB  |  138 lines
  1. function [b,a] = fir2(nn, ff, aa, npt, lap, wind)
  2. %FIR2    FIR filter design using the window method - arbitrary filter shape.
  3. %    B = FIR2(N,F,M) designs an N'th order FIR digital filter with the
  4. %    frequency response specified by vectors F and M, and returns the
  5. %    filter coefficients in length N+1 vector B.  Vectors F and M specify
  6. %    the frequency and magnitude breakpoints for the filter such that
  7. %    PLOT(F,M) 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. %    By default FIR2 uses a Hamming window.  Other available windows,
  13. %    including Boxcar, Hanning, Bartlett, Blackman, Kaiser and Chebwin
  14. %    can be specified with an optional trailing argument.  For example,
  15. %    B = FIR2(N,F,M,bartlett(N+1)) uses a Bartlett window.
  16. %    B = FIR2(N,F,M,chebwin(N+1,R)) uses a Chebyshev window.
  17. %
  18. %    See also FIR1, BUTTER, CHEBY1, CHEBY2, YULEWALK, FREQZ and FILTER.
  19.  
  20. %    Author: L. Shure 3-27-87
  21. % Revised: C. Denham 7-26-90
  22. %    (c) Copyright 1987, by The MathWorks, Inc.
  23.  
  24. %    Optional input arguments (see the user's guide):
  25. %      npt - number for interpolation
  26. %      lap - the number of points for jumps in amplitude
  27.  
  28. nn = nn + 1;
  29. if (nargin < 3 | nargin > 6)
  30.    error('Wrong number of input parameters')
  31. end
  32. if (nargin > 3)
  33.     if nargin == 4
  34.         if length(npt) == 1
  35.            if (2 .^ round(log(npt)/log(2))) ~= npt
  36.             % NPT is not an even power of two
  37.               npt = 2.^ceil(log(npt)/log(2));
  38.            end 
  39.         else
  40.             wind = npt;
  41.             npt = 512;
  42.         end
  43.         lap = fix(npt/25);
  44. % wind = hamming(nn);
  45.     elseif nargin == 5
  46.         if length(npt) == 1
  47.             if (2 .^ round(log(npt)/log(2))) ~= npt
  48.             % NPT is not an even power of two
  49.                 npt = 2.^ceil(log(npt)/log(2));
  50.             end 
  51.             if length(lap) == 1
  52.                 wind = hamming(nn);
  53.             else
  54.                 wind = lap;
  55.                 lap = fix(npt/25);
  56.             end
  57.         else
  58.             wind = npt;
  59.             npt = lap;
  60.             lap = fix(npt/25);
  61.         end
  62.     end
  63. elseif nargin == 3
  64.     if nn < 1024
  65.         npt = 512;
  66.     else
  67.         npt = 2.^ceil(log(nn)/log(2));
  68.     end
  69.     wind = hamming(nn);
  70.     lap = fix(npt/25);
  71. end
  72.  
  73. if nn ~= length(wind)
  74.     error('The specified window must be the same as the filter length')
  75. end
  76.    
  77. [mf,nf] = size(ff);
  78. [ma,na] = size(aa);
  79. if ma ~= mf | na ~= nf
  80.    error('You must specify the same number of frequencies and amplitudes')
  81. end
  82. nbrk = max(mf,nf);
  83. if mf < nf
  84.    ff = ff';
  85.    aa = aa';
  86. end
  87.  
  88. if abs(ff(1)) > eps | abs(ff(nbrk) - 1) > eps
  89.    error('The first frequency must be 0 and the last 1')
  90. end
  91.  
  92. % interpolate breakpoints onto large grid
  93.   
  94. H = zeros(1,npt);
  95. nint=nbrk-1;
  96. df = diff(ff');
  97. if (any(df < 0))
  98.    error('Frequencies must be non-decreasing')
  99. end
  100.  
  101. npt = npt + 1;   % Length of [dc 1 2 ... nyquist] frequencies.
  102.  
  103. nb = 1;
  104. H(1)=aa(1);
  105. for i=1:nint
  106.     if df(i) == 0
  107.        nb = nb - lap/2;
  108.        ne = nb + lap;
  109.     else
  110.        ne = fix(ff(i+1)*npt);
  111.     end
  112.     if (nb < 0 | ne > npt)
  113.        error('Too abrupt an amplitude change near end of frequency interval')
  114.     end
  115.     j=nb:ne;
  116.     if nb == ne
  117.         inc = 0;
  118.     else
  119.         inc = (j-nb)/(ne-nb);
  120.     end
  121.     H(nb:ne) = inc*aa(i+1) + (1 - inc)*aa(i);
  122.     nb = ne + 1;
  123. end
  124.  
  125. % Fourier time-shift.
  126.  
  127. dt = 0.5 .* (nn - 1);
  128. rad = -dt .* sqrt(-1) .* pi .* (0:npt-1) ./ (npt-1);
  129. H = H .* exp(rad);
  130.  
  131. H = [H conj(H(npt-1:-1:2))];   % Fourier transform of real series.
  132. ht = real(ifft(H));            % Symmetric real series.
  133.  
  134. b = ht(1:nn);         % Raw numerator.
  135. b = b .* wind(:).';   % Apply window.
  136.  
  137. a = 1;                % Denominator.
  138.