fir.m

matlab算法集 matlab算法集

function b = fir (g,p)
%-----------------------------------------------------------------------
% Usage:       b = fir (g,p)
%
% Description: Design an FIR linear-phase filter with speciefied
%              magnitude response.  Use the Fourier series method 
%              with a Hamming window to reduce the side lobes.
%
% Inputs:      g = string containing name of function which specifies
%                  the desired magnitude response of the filter.
%                  The function g should be of the form
%
%                       function a = g(f)
%
%                  The scalar input 0 <= f <= 1 represents the 
%                  normalized frequency.  When g is called with
%                  input f, it must return the desired filter 
%                  gain a = g(f).
%              p = filter order (p >= 0)
% 
% Outputs:     b = (2p + 1) by 1 filter coefficient vector.  The
%                  FIR filter output is:
%
%                  y(k) = b[1]*u[k] + b[2]*u[k-1] + ... + b[2p+1]*u[k-2p]  
%   
% Note:        The normalized frequency f is the fraction of the 
%              Nyquist frequency fn = fs/2 where fs is the sampling
%              frequency.
%-----------------------------------------------------------------------
   
% Initialize

   p = args (p,0,p,2,'fir');
   b = zeros (2*p+1,1);
   c = zeros (p,1);
   n = 50;
   h = 1/n;

% Compute filter coefficients using Simpson integration

   disp ('Computing filter coefficients ... ')
   b(p+1) = simpson (0,1,n,g);
   for k = 1 : p
      c(k) = feval(g,0) + cos(pi*k)*feval(g,1);
      for i = 1 : n-1
         if mod(i,2)
            c(k) = c(k) + 4*cos(pi*k*(i*h))*feval(g,i*h);
         else
            c(k) = c(k) + 2*cos(pi*k*(i*h))*feval(g,i*h); 
         end
      end  
      c(k) = h*c(k)/3;
   end

% Add Hamming window

   for k = 1 : p
      b(p+1+k) = (0.54 + 0.46*cos(k*pi/p))*c(k); 
      b(p+1-k) = b(p+1+k);
   end
%-----------------------------------------------------------------------