/scilab-master-1333395999/modules/signal_processing/macros/pspect.sci

// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - 1988 - C. Bunks
// 
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at    
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function [sm,cwp]=pspect(sec_step,sec_leng,wtype,x,y,wpar)
//[sm [,cwp] ]=pspect(sec_step,sec_leng,wtype,x [,y] [,wpar])
//Cross-spectral estimate between x and y if both are given
//and auto-spectral estimate of x otherwise.
//Spectral estimate obtained using the modified periodogram method
// x        :
//            - if vector, the first signal data 
//            - amount of input data if scalar. In this case the first signal
//              is read by batch by a   user defined function xbatch=getx(n,offset).)
// y        : optional, if present cross correlation is made else auto correlation
//            - if vector, the second  signal data .It must have the same number of element than x
//            - if scalar, the second signal
//              is read by batch by a  user defined function ybatch=gety(n,offset).)
// sec_step : offset of each data window. 
//            if sec_step==sec_leng/2 50% overlap is made
// sec_leng : Number of points of the window 
// wtype    : window type ('re','tr','hm','hn','kr','ch')
// wpar     : optional window parameters
//          :         for wtype='kr', wpar>0
//          :         for wtype='ch', 0<wpar(1)<.5, wpar(2)>0
// sm       : power spectral estimate in the interval [0,1]
// cwp      : unspecified Chebyshev window parameter
//!

  [lhs,rhs]=argn(0);
  
  
  if sec_step>=sec_leng then
    error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Argument #%d expected to be less than argument #%d.\n"),"pspect",1,2,1,2));
  end
  if and(wtype<>['re','tr','hm','hn','kr','ch']) then
    error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
                   "pspect",3,"''re'',''tr'',''hm'',''hn'',''kr'',''ch''"));
  end
  
  //Analyze calling sequence and construct window
  if rhs==4 then, 
    //pspect(sec_step,sec_leng,wtype,x)
    w=window(wtype,sec_leng);
    cross=%f; //autocorrelation
  elseif rhs==5 then,
    //pspect(sec_step,sec_leng,wtype,x,y) or pspect(sec_step,sec_leng,wtype,x,wpar)
    if wtype=='kr' then,//pspect(sec_step,sec_leng,'kr',x,wpar)
      wpar=y;
      w=window(wtype,sec_leng,wpar);cwp=[];
      cross=%f; //autocorrelation
    elseif wtype=='ch' then,//pspect(sec_step,sec_leng,'ch',x,wpar)
      wpar=y;
      [w,cwp]=window(wtype,sec_leng,wpar);
      cross=%f; //autocorrelation
    else,//pspect(sec_step,sec_leng,wtype,x,y)
      w=window(wtype,sec_leng);cwp=[];
      cross=%t; //cross-correlation
    end,
  else //pspect(sec_step,sec_leng,wtype,x,y,wpar) 
    [w,cwp]=window(wtype,sec_leng,wpar);
    cross=%t; //cross-correlation
  end,
  wpower=w*w';//the window energy
  
  //Make x and y row vectors
  x=matrix(x,1,-1);
  if cross then 
    y=matrix(y,1,-1);
    if size(x,'*')<>size(y,'*') then
      error(msprintf(gettext("%s: Arguments #%d and #%d must have the same sizes.\n"),"pspect",4,5));
    end
  end
  
  //Average periodograms
  sm=0*w;
  if size(x,'*')==1 then  //Batch processing of x and y data 
    
    //get number of sections to be used                 
    nsecs=int((x-sec_leng+sec_step)/sec_step);//x  contains the number of points of the signals 
    ovrlp=sec_leng-sec_step;
    xd=[0*ones(1,sec_step) getx(ovrlp,1)];
    if cross then
      yd=[0*ones(1,sec_step) gety(ovrlp,1)];
    end
    for k=1:nsecs
      xd(1:ovrlp)=xd(sec_step+1:sec_leng);
      xd(ovrlp+1:sec_leng)=getx(sec_step,sec_leng+(k-1)*sec_step+1);
      xw=w.*(xd-(sum(xd)/sec_leng));
      fx=fft(xw,-1);
      if cross then
        yd(1:ovrlp)=yd(sec_step+1:sec_leng);
        yd(ovrlp+1:sec_leng)=gety(sec_step,sec_leng+(k-1)*sec_step+1);
        yw=w.*(yd-(sum(yd)/sec_leng));
        sm=sm+fx.*conj(fft(yw,-1));
      else
        sm=sm+real(fx.*conj(fx));
      end
    end
  else // Signal data given in x and y if cross-correlation
    
    //get number of sections to be used
    nsecs=int((size(x,'*')-sec_leng+sec_step)/sec_step);
    ind=(1:sec_leng);
    for k=1:nsecs
      indd=ind+(k-1)*sec_step*ones(1:sec_leng);
      xd=x(indd);
      xe=w.*(xd-(sum(xd)/sec_leng));
      fx=fft(xe,-1);
      if cross then
        yd=y(indd);
        ye=w.*(yd-(sum(yd)/sec_leng));
        sm=sm+fx.*conj(fft(ye,-1));
      else
        sm=sm+real(fx.*conj(fx));
      end
    end
  end
  
  //Normalization
  sm=sm/(nsecs*wpower);
endfunction