firkaiser.m

Matlab时频分析工具箱,希望能对大家有所帮助啊

function g=firkaiser(p1,p2,p3,p4)
%FIRKAISER  Kaiser-Bessel window
%   Usage:  g=firkaiser(L,beta);
%           g=firkaiser(L,beta,centering);
%
%   FIRKAISER(L,beta) computes the Kaiser-Bessel window of length L with
%   parameter beta. The smallest element of the window is set to zero when
%   the window has an even length. This gives the window perfect whole-point
%   even symmetry, and makes it possible to use the window for a Wilson
%   basis.
%
%   FIRKAISER(L,beta,centering) will create a window centered as 
%   specified by centering. The default (centering=0) is to return
%   a whole-point centered window, while a value of .5 will produce a
%   half-point even function (traditional signal processing style).
%
%   SEE ALSO: FIRWIN
%
%   REFERENCES:
%     A. V. Oppenheim and R. W. Schafer. Discrete-time signal processing.
%     Prentice Hall, Englewood Cliffs, NJ, 1989.
%     

% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

error(nargchk(2,4,nargin));

if ischar(p1)
  % The user supplied a type argument

  if nargin==2
    error('Two few input arguments.');
  end;

  stype=p1;
  L=p2;
  beta=p3;
  
  if nargin==3
    centering=0;
  else
    centering=p4;
  end;

else

  stype='normal';
  L=p1;
  beta=p2;
  
  if nargin==2
    centering=0;
  else
    centering=p3;
  end;

end;

if numel(beta)>1
  error('beta must be a scalar.');
end;

m=2*((0:L-1)+centering).'/L-1;
switch lower(stype)
  case {'n','normal'}
    g=besseli(0,beta*sqrt(1-m.^2))/besseli(0,beta);
    g=fftshift(g);

    if rem(L,2)==0
      % Explicitly zero last element.
      g(L/2+1)=0;
    end;
    
  case {'d','deriv','derived'}
    g=besseli(0,beta*sqrt(1-m.^2))/besseli(0,beta);    
    g1=sqrt(cumsum(g(1:L/2))./sum(g(1:L/2)));

    if centering==0
      g=[flipud(g1);...
	 0;
	 g1(1:L/2-1)];
    else
      g=[flipud(g1);...
	 g1];
    end;
          
  otherwise
    error('Unknown window type.');
end;