Spectral_Analysis :
    This function performs a spatial spectral analysis with different options on a time series of SLA profiles.
    
    :parameter dx: sampling distance
    :parameter Ain: 2D table of sla data with time along 2nd axis (NXxNT with NX the spatial length and NT the time length)
    :keyword tapering: apply tapering to the data
    
      * If this keyword is of type bool : apply hamming window.
      * If this keyword is a string : apply a hamming ('hamm'), hann ('hann'), kaiser-bessel ('kaiser'), kaiser-bessel ('blackman') or no ('none') tapering function.
      * If this keyword is an :class:`numpy.array` object : apply this array as taper.
       
    :keyword overlap: overlap coefficient of the windows (0.75 means 75% overlap).
    :keyword wsize: size of the sub-segments.
    :keyword normalise: If True, normalise the spectrum by its overall energy content.
    :keyword detrend: If True, removes a linear trend to the segmented signal (if tapered) or to the whole signal (if not tapered).
    :keyword integration: If True, integrate the spectrum between 2 frequencies.
    :keyword alpha: used to compute the input (beta) of the kaiser-bessel taper.
    :keyword ARspec: Applies an Auto-Regression model of the order provided as value of this parameter. 
    
    :return: a spectrum structure
    
       .. code-block:: python
       
          {'esd':esd,       #Energy Spectral Density
           'psd':psd,       #Power Spectral Density
           'fq':fq,         #frequency
           'p':p,           #wavelength
           'params':params} #tapering parameters.
    
    :author: Renaud DUSSURGET (RD) - LER/PAC, Ifremer
    :change: Created by RD, December 2012

def spectral_analysis(dx,Ain,tapering=None,overlap=None,wsize=None,alpha=3.0,detrend=False,normalise=False,integration=True,average=True,ARspec=None):
    """
    Spectral_Analysis :
    This function performs a spatial spectral analysis with different options on a time series of SLA profiles.
    
    :parameter dx: sampling distance
    :parameter Ain: 2D table of sla data with time along 2nd axis (NXxNT with NX the spatial length and NT the time length)
    :keyword tapering: apply tapering to the data
    
      * If this keyword is of type bool : apply hamming window.
      * If this keyword is a string : apply a hamming ('hamm'), hann ('hann'), kaiser-bessel ('kaiser'), kaiser-bessel ('blackman') or no ('none') tapering function.
      * If this keyword is an :class:`numpy.array` object : apply this array as taper.
       
    :keyword overlap: overlap coefficient of the windows (0.75 means 75% overlap).
    :keyword wsize: size of the sub-segments.
    :keyword normalise: If True, normalise the spectrum by its overall energy content.
    :keyword detrend: If True, removes a linear trend to the segmented signal (if tapered) or to the whole signal (if not tapered).
    :keyword integration: If True, integrate the spectrum between 2 frequencies.
    :keyword alpha: used to compute the input (beta) of the kaiser-bessel taper.
    :keyword ARspec: Applies an Auto-Regression model of the order provided as value of this parameter. 
    
    :return: a spectrum structure
    
       .. code-block:: python
       
          {'esd':esd,       #Energy Spectral Density
           'psd':psd,       #Power Spectral Density
           'fq':fq,         #frequency
           'p':p,           #wavelength
           'params':params} #tapering parameters.
    
    :author: Renaud DUSSURGET (RD) - LER/PAC, Ifremer
    :change: Created by RD, December 2012
    """
    
    A=Ain.copy()

    #Check dimensions
    sh = A.shape
    ndims = len(sh)
    N = sh[0] #Time series are found along the last dimension
    
    #If vector, add one dimension
    if ndims == 1 :
        A = A.reshape((N,1))
        sh = A.shape
        ndims = len(sh)
    
    nr = sh[1] #Number of repeats  
    nt = nr
    
#    gain=1.0 #Scaling gain... (used for tapering issues)
    
#    #Get the overall energy
#    spec=get_spec(dx, A[:,0])
#    F=spec['fq']   
#    Eref = ((A[:,0]-A[:,0].mean())**2).sum() #Get the reference energy level
#    ESDref=spec['esd']
#    SFactor=Eref/spec['esd'].sum()
#    ESDref*=SFactor
#    PSDref=spec['psd']*SFactor
#    print 'Check parseval theorem : SUM|Y(f)|²={0}, SUM|y(t)|²={1}'.format(spec['esd'].sum(),((A[:,0]-A[:,0].mean())**2).sum())
    
    #Apply tapering if asked
    ########################
    if tapering is not None:
        
        #Set tapering defaults
        overlap=0 if overlap is None else overlap
        wsize=N if wsize is None else wsize

        #Get time splitting (tapering) parameters
        #########################################
        a = np.float32(wsize)
        b = np.float32(overlap) 
        c = np.float32(N) 
        nn=np.floor((c - (a * b))/(a - (a * b))) #This is the number of segments
        print 'Number of windows :{0}\nTotal windowed points : {1} ({2} missing)\nTotal points : {3}'.format(nn,nn*wsize,N - nn*wsize*overlap,N)
        
        ix = np.arange(nn) * ((1.0 - b) * a) #These are the starting points of each segments

        #Moving window
        ##############
        dum = np.zeros((wsize, nn, nr),dtype=np.float64)
        for j in np.arange(nr):
            for i in np.arange(nn): #looping through time to get splitted time series 
                dum[:,i,j] = detrend_fun(np.arange(wsize),A[ix[i] : ix[i] + wsize,j]) if detrend else A[ix[i] : ix[i] + wsize,j]
        
        #Set up tapering window
        #######################
        beta=np.pi*alpha
        hamm = np.hamming(wsize)
        hann = np.hanning(wsize)
        kbess = np.kaiser(wsize,beta)
        blackman = np.blackman(wsize)
        notaper = np.ones(wsize) #overpass tapering option
        gain=1.0
        
        
        if isinstance(tapering,bool) : which='hamm'
        elif isinstance(tapering,str) :
            if tapering.upper() == 'HAMMING' :
                which='hamm'
                gain=np.sum(hamm)/wsize #0.530416666667
            elif tapering.upper() == 'HANNING' :
                which='hann'
                gain=np.sum(hann)/wsize #0.489583333333
            elif tapering.upper() == 'KAISER' :
                which='kbess'
                gain=np.sum(kbess)/wsize #0.394170357504
            elif tapering.upper() == 'NONE' :
                which='notaper'
                gain=1.0
            elif tapering.upper() == 'BLACKMAN' :
                which='blackman'
                gain=np.sum(blackman)/wsize
            else : raise Exception('Unknown taper {0}'.format(tapering))
        elif isinstance(tapering,np.ndarray) :
            which='tapering'
            gain=np.sum(tapering)/wsize
        else :
            raise Exception('Bad value for tapering keyword')
        exec('window='+which)
        window = np.repeat(window,nn*nr).reshape((wsize,nn,nr))
    
        #Apply tapering on segmented data
        A=dum.copy()*window
        A=A.reshape(wsize,nr*nn) #Reshapa matrix
        nr=nn*nr
    else :
        if detrend :
            for i in np.arange(nr): A[:,i] = detrend_fun(np.arange(N),A[:,i]) if detrend else A[:,i]
        gain=1.0
    
    #Run transform
    ###############
    for i in np.arange(nr):
        if ARspec is not None : spec=yule_walker_regression(dx, A[:,i], ARspec)
        else : spec=get_spec(dx, A[:,i],integration=integration,gain=gain)
        if i == 0:
            esd = spec['esd']
            psd = spec['psd']
            fq = spec['fq']
            phase = spec['phase']
            if ARspec: model =spec['psd'].model
        else : 
            esd = np.append(esd,spec['esd'])
            psd = np.append(psd,spec['psd'])
            phase = spec['phase']
            if ARspec: model =np.append(model,spec['psd'].model)
    
#    factor=((A[:,0]-A[:,0].mean())**2).sum()/spec['esd'].sum()
    #Average spectrum
    #################
    nf=len(fq)
    p=1./fq
    esd=esd.reshape(nr,nf)
    psd=psd.reshape(nr,nf)
    
    if average :
        esd=(np.sum(esd,axis=0)/nr)#/gain
        psd=(np.sum(psd,axis=0)/nr)#/gain
    
    psd = psd * (gain**0.5)
#    print gain, np.sqrt(gain), gain **2, gain*0.5, gain/2.
#    esd=(np.sum(esd,axis=0))#/gain
#    psd=(np.sum(psd,axis=0))#/gain

    if ARspec:
        mparam,msig=zip(*tuple([ (el['parameters'],el['sig']) for el in model]))
        deg=ARspec
        mparam=np.concatenate(mparam).reshape((len(model),deg))
        msig=np.array(msig)
        
        if average :
            mparam=mparam.mean(axis=0)
            msig=np.nansum(msig)/msig.size
            
        mstr={'parameters':mparam,'sig':msig,'n':N,'deg':ARspec}
        setattr(esd,'model',mstr)
        setattr(psd,'model',mstr)
        
        
    #Normalise by energy content    
    Scaling_Factor=len(fq)/esd.sum()
    if normalise :
        esd*=Scaling_Factor
        psd*=Scaling_Factor
    
    if tapering is not None : return {'params':{'tapering':tapering is not None,'which':which,'wsize':int(wsize),'nwind':int(nn),'overlap':int(100.*overlap),'gain':gain},'psd':psd,'esd':esd,'fq':fq,'p':p,'phase':phase}
    else : return {'params':{'tapering':tapering is not None},'psd':psd,'esd':esd,'fq':fq,'p':p,'phase':phase}