def THDN(signal, sample_rate, filename):

    signal = signal - mean(signal)  # this is so that the DC offset is not the peak in the case of PDM
    windowed = signal * blackmanharris(len(signal)) 

    

    # Find the peak of the frequency spectrum (fundamental frequency), and filter 
    # the signal by throwing away values between the nearest local minima
    f = rfft(windowed)

    #limit the bandwidth
    if len(f) > 24000:
       f[24000:len(f)] = 0

    bandwidth_limited = irfft(f)
    total_rms = rms_flat(bandwidth_limited)

    #for i in range(6000):
       #print abs(f[i])
    i = argmax(abs(f))

    #This will be exact if the frequency under test falls into the middle of an fft bin 
    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed)))
    lowermin, uppermin = find_range(abs(f), i)


    #notch out the signal
    f[lowermin: uppermin] = 0



    noise = irfft(f)
    THDN = rms_flat(noise) / total_rms
    print "THD+N:     %.12f%% or %.12f dB" % (THDN * 100, 20 * log10(THDN))

def filtered_cross_corr(signal1,signal2,bins,smoothing=10,timestep=1):
    """Get the cross-correlation between the signals, first filtering
    the fourier transforms (smoothed top-hat), chopping the fourier
    signals into "bins" bins"""
    signal1 -= signal1.mean()
    signal2 -= signal2.mean()
    x1 = rfft(signal1)
    x2 = rfft(signal2)
    assert len(x1)==len(x2)
    startfreq = arange(1,len(x1),len(x1)/(bins+1))
    position = arange(len(x1))
    freq = fftfreq(len(x1),timestep)
    freqout = zeros(bins)*1.
    out = zeros((bins,len(signal1)))*1.0
    att = ones(len(x1))*1.
    for i in range(bins):
        att[:startfreq[i]] = 0
        att[startfreq[i]:startfreq[i+1]] = 1
        att[startfreq[i+1]:] = 0
        freqout[i] = mean(freq*att[:len(freq)])
        att = smooth(att,smoothing)
        att[0] = 0
        x1dash = x1*att
        sig1dash = irfft(x1dash,len(signal1))
        x2dash = x2*att
        sig2dash = irfft(x2dash,len(signal2))
        out[i] = correlate(sig1dash,sig2dash,'same')
    lag = arange(-len(x2),len(x2),1.)*timestep
    return lag,freqout,out        

    def FFT(self):
        for o in self.object_list:
            x_list = []
            y_list = []
            # print o.chaincode
            print "----"
            print o.object_name + "_1 " + str(o.point_list[0][0]) + " " + str(o.point_list[0][1]) + " 1 1"
            for d in o.chaincode:
                print d,
            print "-1"
            print "----"
            x = y = 0
            for d in o.chaincode:
                x += DIRECTION_MATRIX[d][0]
                y -= DIRECTION_MATRIX[d][1]
                x_list.append(x)
                y_list.append(y)
            n_point = len(x_list)
            x_fft_result = fft.rfft(x_list)
            y_fft_result = fft.rfft(y_list)
            for i in range(20):
                x = x_fft_result[i]
                y = y_fft_result[i]
                print "%e %e %e %e" % ( x.real, x.imag, y.real, y.imag )

            x_list_2 = fft.irfft(x_fft_result[0:10], n_point)
            y_list_2 = fft.irfft(y_fft_result[0:10], n_point)

            print "x"
            print x_list
            print x_list_2
            print "y"
            print y_list
            print y_list_2

def filter_xyz(data): 
	xdata=[a[0] for a in data]
	ydata=[a[1] for a in data]
	zdata=[a[2] for a in data]
	samples=len(ydata)
	xdata=array(xdata)
	ydata=array(ydata)
	zdata=array(zdata)
	try:
#----------------------------------
		xfft=(fft.rfft(xdata))
		yfft=(fft.rfft(ydata))
		zfft=(fft.rfft(zdata))
	

#-------------filtering part --------------
		cutoff=samples/3*2
		xfft=ffilter(xfft,cutoff)	
		yfft=ffilter(yfft,cutoff)
		zfft=ffilter(zfft,cutoff)

		nxdata=fft.irfft(xfft)	
		nydata=fft.irfft(yfft)
		nzdata=fft.irfft(zfft)
	except:
		raise ValueError('null value')
	size=len(nxdata)
	data=[[nxdata[i],nydata[i],nzdata[i]] for i in range(size)]
	return data

 def to_time(self, data_type):
     parity = self.pulses_length % 2
     delta_t = 0.5 / ((self.pulses_length - parity) * self.delta_f)
     data1_time = irfft(self._data[self._data1])
     data2_time = irfft(self._data[self._data2])
     time_pulses_length = data1_time.shape[1]
     pulses_time = data_type(time_pulses_length, self.pulses_nb, delta_t)
     pulses_time._data['Valim'] = self.valim
     pulses_time._data[self._data1] = data1_time
     pulses_time._data[self._data2] = data2_time
     return pulses_time

    def single_step_propagation(self):
        """
        Perform single step propagation. The final Wigner functions are not normalized.
        """
        ################ p x -> theta x ################
        self.wigner_ge = fftpack.fft(self.wigner_ge, axis=0, overwrite_x=True)
        self.wigner_g = fftpack.fft(self.wigner_g, axis=0, overwrite_x=True)
        self.wigner_e = fftpack.fft(self.wigner_e, axis=0, overwrite_x=True)

        # Construct T matricies
        TgL, TgeL, TeL = self.get_T_left(self.t)
        TgR, TgeR, TeR = self.get_T_right(self.t)

        # Save previous version of the Wigner function
        Wg, Wge, We = self.wigner_g, self.wigner_ge, self.wigner_e

        # First update the complex valued off diagonal wigner function
        self.wigner_ge = (TgL*Wg + TgeL*Wge.conj())*TgeR + (TgL*Wge + TgeL*We)*TeR

        # Slice arrays to employ the symmetry (savings in speed)
        TgL, TgeL, TeL = self.theta_slice(TgL, TgeL, TeL)
        TgR, TgeR, TeR = self.theta_slice(TgR, TgeR, TeR)
        Wg, Wge, We = self.theta_slice(Wg, Wge, We)

        # Calculate the remaning real valued Wigner functions
        self.wigner_g = (TgL*Wg + TgeL*Wge.conj())*TgR + (TgL*Wge + TgeL*We)*TgeR
        self.wigner_e = (TgeL*Wg + TeL*Wge.conj())*TgeR + (TgeL*Wge + TeL*We)*TeR

        ################ Apply the phase factor ################
        self.wigner_ge *= self.expV
        self.wigner_g *= self.expV[:(1 + self.P_gridDIM//2), :]
        self.wigner_e *= self.expV[:(1 + self.P_gridDIM//2), :]

        ################ theta x -> p x ################
        self.wigner_ge = fftpack.ifft(self.wigner_ge, axis=0, overwrite_x=True)
        self.wigner_g = fft.irfft(self.wigner_g, axis=0)
        self.wigner_e = fft.irfft(self.wigner_e, axis=0)

        ################ p x  ->  p lambda ################
        self.wigner_ge = fftpack.fft(self.wigner_ge, axis=1, overwrite_x=True)
        self.wigner_g = fft.rfft(self.wigner_g, axis=1)
        self.wigner_e = fft.rfft(self.wigner_e, axis=1)

        ################ Apply the phase factor ################
        self.wigner_ge *= self.expK
        self.wigner_g *= self.expK[:, :(1 + self.X_gridDIM//2)]
        self.wigner_e *= self.expK[:, :(1 + self.X_gridDIM//2)]

        ################ p lambda  ->  p x ################
        self.wigner_ge = fftpack.ifft(self.wigner_ge, axis=1, overwrite_x=True)
        self.wigner_g = fft.irfft(self.wigner_g, axis=1)
        self.wigner_e = fft.irfft(self.wigner_e, axis=1)

    def output(self):
        """ """

        # One dimension
        if len(self._source.shape) == 1:
            source = self._actual_source
            # Use FFT convolution
            if self._fft:
                if not self._toric:
                    P = rfft(source,self._fft_shape[0])*self._fft_weights
                    R = irfft(P, self._fft_shape[0]).real
                    R = R[self._fft_indices]
                else:
                    P = rfft(source)*self._fft_weights
                    R = irfft(P,source.shape[0]).real

                # if self._toric:
                #     R  = ifft(fft(source)*self._fft_weights).real
                # else:
                #     n = source.shape[0]
                #     self._src_holder[n//2:n//2+n] = source
                #     R = ifft(fft(self._src_holder)*self._fft_weights)
                #     R = R.real[n//2:n//2+n]
            # Use regular convolution
            else:
                R = convolve1d(source, self._weights[::-1], self._toric)
            if self._src_rows is not None:
                R = R[self._src_rows]
            return R.reshape(self._target.shape)
        # Two dimensions
        else:
            source = self._actual_source
            # Use FFT convolution
            if self._fft:
                if not self._toric:
                    P = rfft2(source,self._fft_shape)*self._fft_weights
                    R = irfft2(P, self._fft_shape).real
                    R = R[self._fft_indices]
                else:
                    P = rfft2(source)*self._fft_weights
                    R = irfft2(P,source.shape).real

            # Use SVD convolution
            else:
                R = convolve2d(source, self._weights, self._USV, self._toric)
            if self._src_rows is not None and self._src_cols is not None:
                R = R[self._src_rows, self._src_cols]
        return R.reshape(self._target.shape)

    def shift(self,p=None,pdot=None,dm=None):
        if dm is None:
            dm = self.current_dm
        if p is None:
            p = self.current_p
        if pdot is None:
            pdot = self.current_pdot

        f = rfft(self.profs,axis=-1)

        # the sample at phase p is moved to phase p-dmdelays
        dmdelays = (psr_utils.delay_from_DM(dm,self.pfd.subfreqs) -
                    psr_utils.delay_from_DM(self.current_dm,self.pfd.subfreqs))/self.original_fold_p
        
        start_times = self.pfd.start_secs
        pdelays = (start_times/self.original_fold_p) * (p-self.current_p)/self.original_fold_p
        pdotdelays = ((pdot-self.current_pdot)*start_times**2/2.)/self.original_fold_p
        
        f *= np.exp((2.j*np.pi)*
            dmdelays[np.newaxis,:,np.newaxis]*
            np.arange(f.shape[-1])[np.newaxis,np.newaxis,:])
        f *= np.exp((2.j*np.pi)*
            (pdelays+pdotdelays)[:,np.newaxis,np.newaxis]*
            np.arange(f.shape[-1])[np.newaxis,np.newaxis,:])

        self.profs = irfft(f)
        self.current_p = p
        self.current_dm = dm
        self.current_pdt = pdot

def deriv(var, periodic=False):
    """Take derivative of 1D array"""
    n = var.size
    if periodic:
        # Use FFTs to take derivatives
        f = rfft(var)
        f[0] = 0.0 # Zero constant term
        if n % 2 == 0:
            # Even n
            for i in arange(1,n/2):
                f[i] *= 2.0j * pi * float(i)/float(n)
            f[-1] = 0.0 # Nothing from Nyquist frequency
        else:
            # Odd n
            for i in arange(1,(n-1)/2 + 1):
                f[i] *= 2.0j * pi * float(i)/float(n)
        return irfft(f)
    else:
        # Non-periodic function
        result = zeros(n) # Create empty array
        if n > 2:
            for i in arange(1, n-1):
                # 2nd-order central difference in the middle of the domain
                result[i] = 0.5*(var[i+1] - var[i-1])
            # Use left,right-biased stencils on edges (2nd order)
            result[0]   = -1.5*var[0]   + 2.*var[1]   - 0.5*var[2]
            result[n-1] =  1.5*var[n-1] - 2.*var[n-2] + 0.5*var[n-3]
        elif n == 2:
            # Just 1st-order difference for both points
            result[0] = result[1] = var[1] - var[0]
        elif n == 1:
            result[0] = 0.0
        return result

def record2vecs(File):
    # Read the rr data (Note that it is sampled irregularly and
    # consists of the r-times in centiseconds) and return a high
    # frequency spectrogram.
    import cinc2000
    from numpy.fft import rfft, irfft
    data = []
    for line in File:
        data.append(float(line)/100)
    File.close()
    # Now data[i] is an r-time in seconds.
    hrd = cinc2000.R_times2Dev(data)
    # hrd are heart rate deviations sampled at 2 Hz
    pad = np.zeros(((Glength+len(hrd)),),np.float64)
    pad[Gl_2:len(hrd)+Gl_2] = hrd
    # Now pad is hrd with Gl_2 zeros before and after
    N_out = len(hrd)//RDt
    result = np.zeros((N_out,Fw/2),np.float64)
    mags = []
    for k in range(N_out):
        i = int(RDt*k)
        WD = Gw*pad[i:i+Glength] # Multiply data by window fuction
        FT = rfft(WD,n=Fw)
        SP = np.conjugate(FT)*FT  # Periodogram
        temp = rfft(SP,n=Fw//2)
        SP = irfft(temp[0:int(0.1*Fw/2)],n=Fw/2)
        # Low pass filter in freq domain.  Pass below 0.1 Nyquist
        temp = SP.real[:Fw/2]
        mag = math.sqrt(np.dot(temp,temp))
        result[k,:] = temp/mag
        mags.append(math.log(mag))
        # result[k,:] is a unit vector and a smoothed periodogram
    return [result,mags]

    def induced_voltage_generation(self, Beam, length = 'slice_frame'):
        '''
        *Method to calculate the induced voltage from the inverse FFT of the
        impedance times the spectrum (fourier convolution).*
        '''
        if self.recalculation_impedance:
            self.sum_impedances(self.frequency_array)
        
        self.slices.beam_spectrum_generation(self.n_fft_sampling)

        
        if self.save_individual_voltages:
            for i in range(self.len_impedance_source_list):
                self.matrix_save_individual_voltages[:,i] = - Beam.charge * e * Beam.ratio * irfft(self.matrix_save_individual_impedances[:,i] * self.slices.beam_spectrum)[0:self.slices.n_slices] * self.slices.beam_spectrum_freq[1] * 2*(len(self.slices.beam_spectrum)-1)
            self.induced_voltage = np.sum(self.matrix_save_individual_voltages,axis=0)     
        else:
            induced_voltage = - Beam.charge * e * Beam.ratio * irfft(self.total_impedance * self.slices.beam_spectrum) * self.slices.beam_spectrum_freq[1] * 2*(len(self.slices.beam_spectrum)-1) 
            self.induced_voltage = induced_voltage[0:self.slices.n_slices]
            
        
        
        if isinstance(length, int):
            max_length = len(induced_voltage)
            if length > max_length:
                induced_voltage = np.lib.pad(induced_voltage, (0, length-max_length), 'constant', constant_values=(0,0))
            return induced_voltage[0:length]

def bandpass_gaussian(data, dt, period, alpha):
    """
    Bandpassing real data (in array *data*) with a Gaussian
    filter centered at *period* whose width is controlled
    by *alpha*:

      exp[-alpha * ((f-f0)/f0)**2],

    with f the frequency and f0 = 1 / *period*.
    *dt* is the sampling interval of the data.

    @type data: L{numpy.ndarray}
    @type dt: float
    @type period: float
    @type alpha: float
    @rtype: L{numpy.ndarray}
    """
    # Fourier transform
    fft_data = rfft(data)
    # aray of frequencies
    freq = rfftfreq(len(data), d=dt)

    # bandpassing data
    f0 = 1.0 / period
    fft_data *= np.exp(-alpha * ((freq - f0) / f0) ** 2)

    # back to time domain
    return irfft(fft_data, n=len(data))

 def autocorr(self, x):
   """ multi-dimensional autocorrelation with FFT """
   X = rfft(x, n=(x.shape[1]*2-1), axis=1)
   xr = irfft(X * X.conjugate(), axis=1).real
   xr = fftshift(xr, axes=1)
   xr = xr.sum(axis=1)
   return xr

    def ifft(self):
        """Compute the one-dimensional discrete inverse Fourier
        transform of this `FrequencySeries`.

        Returns
        -------
        out : :class:`~gwpy.timeseries.TimeSeries`
            the normalised, real-valued `TimeSeries`.

        See Also
        --------
        :mod:`scipy.fftpack` for the definition of the DFT and conventions
        used.

        Notes
        -----
        This method applies the necessary normalisation such that the
        condition holds:

        >>> timeseries = TimeSeries([1.0, 0.0, -1.0, 0.0], sample_rate=1.0)
        >>> timeseries.fft().ifft() == timeseries
        """
        from ..timeseries import TimeSeries
        nout = (self.size - 1) * 2
        # Undo normalization from TimeSeries.fft
        # The DC component does not have the factor of two applied
        # so we account for it here
        dift = npfft.irfft(self.value * nout) / 2
        new = TimeSeries(dift, epoch=self.epoch, channel=self.channel,
                         unit=self.unit, dx=1/self.dx/nout)
        return new

def cross_correlate(histogram, template, n_harmonics=None, upsample=16):
    """Find the shift required to align histogram with template

    """
    n = max(len(histogram),len(template))

    h_ft = rfft(histogram)
    t_ft = rfft(template)

    if len(h_ft)<len(t_ft):
        h_ft = concatenate((h_ft,zeros(len(t_ft)-len(h_ft))))
    elif len(t_ft)<len(h_ft):
        t_ft = concatenate((t_ft,zeros(len(h_ft)-len(t_ft))))

    if n_harmonics is not None:
        h_ft[n_harmonics+1:]*=0.
        t_ft[n_harmonics+1:]*=0.
    h_ft[0] = 0
    t_ft[0] = 0

    cross_correlations = irfft(conjugate(h_ft)*t_ft,n*upsample)
    shift = argmax(cross_correlations)/float(len(cross_correlations))

    assert 0<=shift<1

    #FIXME: warn if double-peaked

    return shift

def filterFromXFer(ds, dpathXFer='/', newpath='/filter', nyquist=5000):
	d = ds.getSubData(dpathXFer)
	fs, start = d.fs(), d.start()
	gain = d.data[:,0]
	phase = d.data[:,1]
	if not nyquist:
		nyquist = start + gain.shape[0]/fs
	if start:
		npad = int(round(start*fs))
		gain = concatenate([zeros(npad), gain])
		phase = concatenate([zeros(npad), phase])
	else:
		gain[0]=0
		phase[0]=0
	c = gain*exp(1j*phase)
	ts = irfft(c, nyquist*2)
	n = int(ts.shape[0]/2.0)
	ts = concatenate([ts[n:], ts[:n]])
	#phase correction by one sample point here. Why is that?
	#it seems that irfft doesn't use an odd number of Fourier points, and so the center of the 
	#spectrum gets shifted by one? Overtly specifying an odd number of transform points doesn't solve 
	#the problem though. using nyquist*2+1 
	ts=ts[1:]
	head = {"SampleType":"timeseries", "SamplesPerSecond":nyquist*2}
	ds.createSubData(newpath, data=ts, head=head, delete=True)
	analyzeFilterWN(ds, dpathFilt=newpath, useWindowedFFT=True)

def STIFT(stft_windows, nfft, nskip):
  "Sum contributions from each window. No normalization."
  r = np.zeros(nskip * (len(stft_windows) - 1) + nfft)
  h = hamming(nfft)
  for w, w_0 in zip(stft_windows, xrange(0, r.size, nskip)):
    r[w_0 : w_0+nfft] += irfft(w, nfft) * h
  return r

def lpc_python(x, order):
    def levinson(R, order):
        """ input: autocorrelation and order, output: LPC coefficients """
        a = np.zeros(order + 2)
        a[0] = -1
        k = np.zeros(order + 1)
        # step 1: initialize prediction error "e" to R[0]
        e = R[0]
        # step 2: iterate over [1:order]
        for i in range(1, order + 1):
            # step 2-1: calculate PARCOR coefficients
            k[i] = (R[i] - np.sum(a[1:i] * R[i - 1 : 0 : -1])) / e
            # step 2-2: update LPCs
            a[i] = np.copy(k[i])
            a_old = np.copy(a)
            for j in range(1, i):
                a[j] -= k[i] * a_old[i - j]
            # step 2-3: update prediction error "e"
            e = e * (1.0 - k[i] ** 2)
        return -1 * a[0 : order + 1], e, -1 * k[1:]

    # N: compute next power of 2
    n = x.shape[0]
    N = int(np.power(2, np.ceil(np.log2(2 * n - 1))))
    # calculate autocorrelation using FFT
    X = rfft(x, N)
    r = irfft(abs(X) ** 2)[: order + 1]
    return levinson(r, order)

  def getview (self, view, pbar):
  # {{{
    import numpy 
    saxis = self.saxis
    # Get bounds of slice on smoothing axis
    ind = view.integer_indices[saxis]
    st, sp = numpy.min(ind), numpy.max(ind)
    # input is the whole range
    insl = slice(0, self.shape[saxis],1)
    # output is the required slice
    outsl = [ ind if i == saxis else slice(None) for i in range(self.naxes)]
    # Get source data
    aview = view.modify_slice(saxis, insl)
    src = aview.get(self.var, pbar=pbar)
    maxharm= self.maxharm
    smsl = [ slice(maxharm,None) if i == saxis else slice(None) for i in range(self.naxes)]
    # calculate harmonics and output required data
    from numpy import fft 
    if 'complex' in self.dtype.name:
      ct=fft.fft(src,self.shape[saxis],saxis)
      smsl=[ slice(maxharm,-maxharm+1) if i == saxis else slice(None) for i in range(self.naxes)] 
      ct[smsl]=0
      st = fft.ifft(ct, self.shape[saxis], saxis)
    else:
      ct=fft.rfft(src,self.shape[saxis],saxis)
      ct[smsl]=0
      st = fft.irfft(ct, self.shape[saxis], saxis)
    


    return st[outsl].astype(self.dtype)

 def macro_mismatch(self,p1,p2):
     """
     Performs double convolution with two different periods to calculate
     macroscopic average of a charge density along the z-axis.
     """
     from numpy.fft import rfft,irfft
     # Convert periods to direct units, if given in cartesian
     if p1 > 1.0:
         p1 = p1/self.unitcell.cell_vec[2,2]
     if p2 > 1.0:
         p2 = p2/self.unitcell.cell_vec[2,2]
     # Create xy-plane averaged density
     micro_z = self.integrate_z_density()
     # Create convolutions
     z_pos = np.linspace(0,1,len(micro_z))
     # Find index of lowest lower bound for p1
     low = 1.-p1/2.
     i1 = len(z_pos)-1
     while True:
         if z_pos[i1] <= low:
             i1 += 1
             break
         i1 -= 1
     #Find index of lowest upper bound for p1
     high = p1/2.
     j1 = 0
     while True:
         if z_pos[j1] >= high:
             j1 -= 1
             break
         j1 += 1
     # Find index of lowest lower bound for p2
     low = 1.-p2/2.
     i2 = len(z_pos)-1
     while True:
         if z_pos[i2] <= low:
             i2 += 1
             break
         i2 -= 1
     #Find index of lowest upper bound for p2
     high = p2/2.
     j2 = 0
     while True:
         if z_pos[j2] >= high:
             j2 -= 1
             break
         j2 += 1
     conv1 = np.zeros(len(micro_z))
     conv1[0:j1+1] = np.ones(j1+1)
     conv1[i1:] = np.ones(len(conv1[i1:]))
     conv2 = np.zeros(len(micro_z))
     conv2[0:j2+1] = np.ones(j2+1)
     conv2[i2:] = np.ones(len(conv2[i2:]))
     # Perform convolutions in Fourier Space
     f_micro_z = rfft(micro_z)
     f_conv1 = rfft(conv1)
     f_conv2 = rfft(conv2)
     f_macro = f_conv2*f_conv1*f_micro_z
     macro_z = irfft(f_macro)/float(np.sum(conv1))/float(np.sum(conv2))
     return macro_z