mywindow.m

基于Matlab的地震数据处理显示和测井数据显示于处理的小程序

function h=mywindow(N,name,param,param2)
%	h=mywindow(n,name,param,param2)
%	yields a window of length N with a given shape.
%
%	n      : length of the window
%	name   : name of the window shape (default : Hamming)
%	param  : optional parameter
%	param2 : second optional parameters
%
%   Possible names are :
%       'Bartlett', 'Blackman', 'Dolph'
%       'Gauss',    'Hamming',  'Hanning',  'Harris', 
%       'none' (same as 'rect') 'Parzen',   'Papoulis', 'rect',
%       'sin'       'triang' (same as 'Bartlett')
%	
%	For the gaussian window, an optionnal parameter "param"
%	sets both end values. The default value is 0.005
%
%   Example : h=mywindow(256,'Gauss',0.005);
%                 figure
%                 plot(h);


if nargin == 0
   error ('At least 1 input parameter is required.')
end

if N <= 0
   error('N should be strictly positive.')
end;

if nargin == 1
   name= 'HAMMING';
else
   name=upper(name);
end

switch name


case 'BARTHANN'
   h=0.38 * (1.0-cos(2.0*pi*(1:N)/(N+1))') ...
   + 0.48 * min((1:N),(N:-1:1))'/(N+1);


case {'BARTLETT','TRIANG'},
   h=2.0*min((1:N),(N:-1:1))'/(N+1);


case 'BLACKMAN'
   ind=(-(N-1)/2:(N-1)/2)' *2.0*pi/N;
   h=0.42 + 0.50*cos(ind) + 0.08*cos(2.0*ind);


case {'DOLPH','DOLF'}
   if rem(N,2) == 0
      oddN=1; 
      N=2*N+1; 
   else 
      oddN=0; 
   end;
   if nargin == 3
      A=10^(param/20); 
   else 
      A=1e-3; 
   end;
   K=N-1; 
   Z0=cosh(acosh(1.0/A)/K); 
   x0=acos(1/Z0)/pi; 
   x=(0:K)/N; 
   indices1=find((x<x0) | (x>1-x0));
   indices2=find((x>=x0) & (x<=1-x0));
   h(indices1)=cosh(K*acosh(Z0*cos(pi*x(indices1))));
   h(indices2)=cos(K*acos(Z0*cos(pi*x(indices2))));
   h=fftshift(real(ifft(A*real(h))));h=h'/h(K/2+1);
   if oddN
      h=h(2:2:K); 
   end;

%{
case 'FEJER'
   nh=fix(N/2);
   ind=(1:nh-1)/nh;
   h1=sin(pi*ind*(nh+1)).^2./(sin(pi*ind).^2*(nh+1));
   h=h1;
   error('Fejer window not yet implemented')
%  h= needs checking
%}

case 'GAUSS'
   if nargin < 3
      param=0.005; 
   end;
   h=exp(log(param) * linspace(-1,1,N)'.^2 );



case 'HAMMING'
   h=0.54 - 0.46*cos(2.0*pi*(1:N)'/(N+1));


case {'HANNING','HANN'}
  h=0.50 - 0.50*cos(2.0*pi*(1:N)'/(N+1));


case 'HARRIS'
   ind=(1:N)' *2.0*pi/(N+1);
   h= 0.35875 ...
     -0.48829 *cos(    ind) ...
     +0.14128 *cos(2.0*ind) ...
     -0.01168 *cos(3.0*ind);


case 'KAISER'
  if nargin == 3
     beta=param; 
  else 
     beta=3.0*pi; 
  end;
  ind=(-(N-1)/2:(N-1)/2)' *2/N;
  h=bessel(0,j*beta*sqrt(1.0-ind.^2))/real(bessel(0,j*beta));


case {'NONE','RECTANG','RECT'} 
   h=ones(N,1);


case 'NUTBESS'
   if nargin == 3
      beta=param; nu=0.5; 
   elseif nargin == 4
      beta=param; 
      nu=param2;
   else 
      beta=3*pi; 
      nu=0.5;
   end;
   ind=(-(N-1)/2:(N-1)/2)' *2/N; 
   h=sqrt(1-ind.^2).^nu .* ...
     real(bessel(nu,j*beta*sqrt(1.0-ind.^2)))/real(bessel(nu,j*beta));


case 'NUTTALL'
 ind=(-(N-1)/2:(N-1)/2)' *2.0*pi/N;
 h= 0.3635819 ...
   +0.4891775*cos(    ind) ...
   +0.1363995*cos(2.0*ind) ...
   +0.0106411*cos(3.0*ind) ;


case 'PAPOULIS'
   ind=(1:N)'*pi/(N+1); 
   h=sin(ind);


case 'PARZEN'
   ind=abs(-(N-1)/2:(N-1)/2)'*2/N; 
   temp=2*(1.0-ind).^3;
   h=min(temp-(1-2.0*ind).^3,temp);

case 'SINE'
   h=sin(pi*(1:N)/(N+1))';

otherwise
   error([' Unknown window: ',name]);

end