GET_SPEC
    :summary: Returns the spectrum of a regularly sampled dataset
    
    :parameter dq: sampling interval (1D)
    :parameter V: data to analyse (1D).
    
    :note: NaN can not be used. 
    
    :return:
    
      * psd: Power Spectral Density
      * esd: Energy Spectral Density
      * fq: frequency
      * p: wavelength (period)
    
    :author: Renaud DUSSURGET (RD) - LER/PAC, Ifremer
    :change: Created by RD, July 2012. Changes
      * 29/08/2012 : Changed the computation of frequencies and the spectral integration (spectrum is averaged at mid-width frequencies)
      * 30/11/2012 : Outstanding changes : corrected the spectral integration for computing psd and corrected the normalisation

def get_spec(dx,Vin,verbose=False,gain=1.0,integration=True):
    """
    GET_SPEC
    :summary: Returns the spectrum of a regularly sampled dataset
    
    :parameter dq: sampling interval (1D)
    :parameter V: data to analyse (1D).
    
    :note: NaN can not be used. 
    
    :return:
    
      * psd: Power Spectral Density
      * esd: Energy Spectral Density
      * fq: frequency
      * p: wavelength (period)
    
    :author: Renaud DUSSURGET (RD) - LER/PAC, Ifremer
    :change: Created by RD, July 2012. Changes
      * 29/08/2012 : Changed the computation of frequencies and the spectral integration (spectrum is averaged at mid-width frequencies)
      * 30/11/2012 : Outstanding changes : corrected the spectral integration for computing psd and corrected the normalisation
    """
    V=Vin.copy()
    N=V.size

    # part 1 - estimate of the wavenumbers
    ######################################
    [k,L,imx] = get_kx(N,dx)
#    print "dx = %s" %(dx)
#    print "kx = ",k[0:3]
#    print "shape",k.shape

    # part 2 - estimate of the spectrum
    # fast fourier transform
    ##########################
    fft=ft.fft(V)/(gain)
    if verbose : print 'Check parseval theorem 1: SUM|Y(f)|²={0}, SUM|y(t)|²={1}'.format(((np.abs(fft)**2)/N).sum(),(V**2).sum()) 

#        Compute magnitude
#        CF.Wilks, 2006. Statistical Methods in the Atmopsheric Sciences. 2nd Edition. International Geophysics Series vol.96. Elsevier.
#        => Ck = [Ak^2 + Bk^2]^(1/2), where e^(iwt) = cos(wt) + i.sin(wt)
    a = fft.real
    b = fft.imag
    c = np.sqrt(a**2.0 + b**2.0) #This is the magnitude
    d = np.arctan(b/a)
    
    if verbose : print 'Check parseval theorem 2: SUM|Y(f)|²={0}, SUM|y(t)|²={1}'.format(((c**2)/N).sum(),(V**2).sum()) 
    
    #Normalise data
    c /= np.float32(N)
    #~ d /= np.float32(N)
    
    if verbose : print 'Check parseval theorem 3: SUM|Y(f)|²={0}, SUM|y(t)|²={1}'.format(((c**2)*N).sum(),(V**2).sum()) 
    
#    mean=fft.real[0]/N #Get average
    
    #Remove negative frequencies
#    phase=np.angle(2*fft[1:imx]) #Get phase (not used yet)
    c = 2*c[1:imx-1] #Multiply by 2 (because of negative frequencies removal) - loses a bit of energy
    d = 2*d[1:imx-1]
    
    if verbose : print 'Check parseval theorem 4: SUM|Y(f)|²={0}, SUM|y(t)|²={1}'.format(((c**2)*(N/2.0)).sum(),((V-V.mean())**2).sum()) 
    
#    if not(N&1 and True or False) : print (np.sum(fft[1:N/2+1])-np.sum(fft[N/2+1:])), fft[N/2] #This is equal to fft[n/2] 
#    cbis=np.absolute(2*fft[0:imx-1]) #Slightly differing from previous one    
#    cbis /= np.float32(N)
   
    # integration of the spectrum
    # shift to have values centered on the intervals
    dk = k[1] - k[0]
    dk_half = dk/2  # half of the interval
        
    k = k[1:imx-1]  
    k_= k[:-1] + dk_half  #Shift wavelengths by half the unit
    
#    cesd = c**2 *(N/2) #Energy (ESD)
    csquared = c ** 2 #Energy (ESD)
    
    if verbose : print 'Check parseval theorem 5: SUM|Y(f)|²={0}, SUM|y(t)|²={1}'.format((csquared*(N/2)).sum(),((V-V.mean())**2).sum())
    
    #Spectral integration
    if integration :
        esd = k_*0.0
        psd = k_*0.0
        for i in np.arange(len(k_)):
            esd[i] = np.sum((csquared * (N/2.0))[(k > (k_[i]-dk)) & (k < (k_[i]+dk))]) / 2.0
            psd[i] = np.sum(csquared[(k > (k_[i]-dk)) & (k < (k_[i]+dk))]) / (2.0**2) #This is variance units integration
        fq=k_
    else :
        esd = csquared
        psd = esd.copy()  /2.
        fq=k.copy()
        
    psd = psd / dk
#    psd = psd/ (N**2.0) / dk    # Normalisation (Danioux 2011) 

    if verbose : print 'Check parseval theorem 6: SUM|Y(f)|²={0}, SUM|y(t)|²={1}'.format(esd.sum(),((V-V.mean())**2).sum())

    #Get frequencies and period  
    p=1/fq
    
    return {'psd':psd,'esd':esd,'fq':fq,'p':p,'gain':gain, 'phase':d}