def Resample(S, ns, nr, p, rate):
  w_s = rfftfreq(ns, 1.0/rate)
  d_s  = 1.*rate / ns
  w_r = rfftfreq(nr, 1.0/rate)
  d_r = 1.*rate / nr
  i_wfloor = (np.floor(w_r[1:]/p / d_s) - 1).astype(np.int)
  i_wceil = (np.ceil(w_r[1:]/p / d_s) - 1).astype(np.int)
  
  w_floor = w_s[1+i_wfloor]
  w_ceil = w_s[1+i_wceil]
  a_floor0, a_ceil0 = (w_r[1:]/p - w_floor), (w_ceil - w_r[1:]/p)
  # Some freqs may be exact matches, e.g. if nr == 2*ns, every other one.
  exact = (a_floor0 == 0.0) & (a_ceil0 == 0.0)
  # Ignore the DC Component and any exact matches in the interpolation.
  R = np.hstack((np.array([S[0]]), np.zeros(w_r.size-1)))
  R[1:][exact] = S[1+i_wfloor[exact]]
  # The rest: interpolate mag/phase from the nearest sampled frequencies.
  absS, angS = np.abs(S[1:]), np.angle(S[1:])
  intrp = ~exact
  # Basic floor/ceil linear interpolation ...
  i_fl, i_cl = i_wfloor[intrp], i_wceil[intrp]
  afl0, acl0 = a_floor0[intrp], a_ceil0[intrp]
  afl = 1.0 - afl0 / (afl0 + acl0)
  acl = 1.0 - acl0 / (afl0 + acl0)
  absR = afl * absS[i_fl] + acl * absS[i_cl]
  angR = afl * angS[i_fl] + acl * angS[i_cl]
  R[1:][intrp] = absR * np.exp(1j * angR)
  assert R.size == w_r.size, (R.size, w_r.size, S.size, ns, nr)
  return R

 def test_definition(self):
     x = [0, 1, 2, 3, 4]
     assert_array_almost_equal(9*fft.rfftfreq(9), x)
     assert_array_almost_equal(9*pi*fft.rfftfreq(9, pi), x)
     x = [0, 1, 2, 3, 4, 5]
     assert_array_almost_equal(10*fft.rfftfreq(10), x)
     assert_array_almost_equal(10*pi*fft.rfftfreq(10, pi), x)

def dfmap_fft(df, indep_var='Time'):
    if df is not None:
        df = df
    new_df = {}
    n = len(df[indep_var])
    a, b = df[indep_var][0], df[indep_var][n-1]
    for key in df:
        if key == indep_var:
            continue
        new_df[key] = rfft(df[key])
        new_df[key] = list(map(np.absolute, new_df[key]))
        new_df[key] = list(map(lambda x: 2*x/n, new_df[key]))

    # avg time between samples
    # assumes uniform sampling
    dt = (b-a) / n
    #dt = (b-a).total_seconds() / n

    # get list of frequencies from the fft
    # new_df['Freq'] = fftfreq(len(self.df[key]), dt)
    new_df['Freq'] = rfftfreq(n, dt)

    #
    new_df = pandas.DataFrame.from_dict(new_df)#.sort_values('Freq')
    #new_df = new_df[abs(new_df.Freq) > 0.4]

    return new_df

    def create_output_signal(self, freq=1000, level=-3, sig_type='sine'):
        '''Creates the array for global out_sig.
        Should eventually handle multiple signal types.
        Currently this only builds sine waves that are periodic
        in our buffer size. It finds the closest frequency to the
        specified frequency.

        takes inputs:
            freq: frequency in Hz.
            level: signal level in dB.
            sig_type: Currently only supports sine waves.
        '''

        if self.out_enable:
            retoggle = True
            self.toggle_out()
        else:
            retoggle = False

        freq_array = fft.rfftfreq(n=self.args.buff_size,
                                  d=(1 / self.args.sample_rate))
        mag_array = np.zeros_like(freq_array)
        closest_freq_index = np.searchsorted(freq_array, freq)
        mag_array[closest_freq_index] = ((10 ** (level / 20)
                                         * self.args.buff_size / 2))
        self.out_sig = np.fft.irfft(a=mag_array, n=self.args.buff_size)
        if retoggle:
            self.toggle_out()

        return (freq_array[closest_freq_index], level, sig_type)

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 decompose(self,l_edges,keep_fourier=False):

		"""
		Decomposes the shear map into its E and B modes components and returns the respective power spectral densities at the specified multipole moments

		:param l_edges: Multipole bin edges
		:type l_edges: array

		:param keep_fourier: If set to True, holds the Fourier transforms of the E and B mode maps into the E and B attributes of the ShearMap instance
		:type keep_fourier: bool. 

		:returns: :returns: tuple -- (l -- array,P_EE,P_BB,P_EB -- arrays) = (multipole moments, EE,BB power spectra and EB cross power)

		>>> test_map = ShearMap.load("shear.fit",format=load_fits_default_shear)
		>>> l_edges = np.arange(300.0,5000.0,200.0)
		>>> l,EE,BB,EB = test_map.decompose(l_edges)

		"""

		#Perform Fourier transforms
		ft_data1 = rfft2(self.data[0])
		ft_data2 = rfft2(self.data[1])

		#Compute frequencies
		lx = rfftfreq(ft_data1.shape[0])
		ly = fftfreq(ft_data1.shape[0])

		#Safety check
		assert len(lx)==ft_data1.shape[1]
		assert len(ly)==ft_data1.shape[0]

		#Compute sines and cosines of rotation angles
		l_squared = lx[np.newaxis,:]**2 + ly[:,np.newaxis]**2
		l_squared[0,0] = 1.0

		sin_2_phi = 2.0 * lx[np.newaxis,:] * ly[:,np.newaxis] / l_squared
		cos_2_phi = (lx[np.newaxis,:]**2 - ly[:,np.newaxis]**2) / l_squared

		#Compute E and B components
		ft_E = cos_2_phi * ft_data1 + sin_2_phi * ft_data2
		ft_B = -1.0 * sin_2_phi * ft_data1 + cos_2_phi * ft_data2

		ft_E[0,0] = 0.0
		ft_B[0,0] = 0.0

		assert ft_E.shape == ft_B.shape
		assert ft_E.shape == ft_data1.shape

		#Compute and return power spectra
		l = 0.5*(l_edges[:-1] + l_edges[1:])
		P_EE = _topology.rfft2_azimuthal(ft_E,ft_E,self.side_angle.to(deg).value,l_edges)
		P_BB = _topology.rfft2_azimuthal(ft_B,ft_B,self.side_angle.to(deg).value,l_edges)
		P_EB = _topology.rfft2_azimuthal(ft_E,ft_B,self.side_angle.to(deg).value,l_edges)

		if keep_fourier:
			self.fourier_E = ft_E
			self.fourier_B = ft_B

		return l,P_EE,P_BB,P_EB

def freqz_fos(b, a, order, nfft, plotfun=None):
    impulse = _create_impulse(nfft)
    response, states = fosfilter(b, a, order, impulse)
    freqresponse = rfft(np.real(response))
    frequencies = rfftfreq(nfft)
    if plotfun:
        plotfun(frequencies, freqresponse)
    return freqresponse, frequencies, response

def fourier_analysis(data, ptree=None):
    time = data['time']
    current = data['current']
    voltage = data['voltage']

    # inspect data
    n = len(time)  # normalization factor for fft
    assert len(current) == n
    assert len(voltage) == n
    d = time[1] - time[0]  # inverse of the sampling rate
    # check sampling spacing is the same everywhere
    for i in range(n - 1):
        assert isclose(time[i + 1] - time[i], d, atol=1e-10, rtol=1e-10)

    # truncate signals
    if ptree:
        steps_per_cycle = ptree.get_int('steps_per_cycle')
        cycles = ptree.get_int('cycles')
        ignore_cycles = ptree.get_int('ignore_cycles')
        assert cycles > ignore_cycles
        assert n == cycles * steps_per_cycle
        time = time[ignore_cycles * steps_per_cycle:]
        current = current[ignore_cycles * steps_per_cycle:]
        voltage = voltage[ignore_cycles * steps_per_cycle:]
    else:
        time = time[int(n / 2):]
        current = current[int(n / 2):]
        voltage = voltage[int(n / 2):]

    n = len(time)
    assert len(current) == n
    assert len(voltage) == n

    if not _is_power_of_two(n):
        warn(
            "(cycles-ignore_cycles)*steps_per_cycles is not a "
            "power of 2 (most efficient for the fourier analysis)",
            RuntimeWarning)

    # perform the actual fourrier analaysis
    fft_current = fft.rfft(current) / n
    fft_voltage = fft.rfft(voltage) / n
    fft_frequency = fft.rfftfreq(n, d)

    # find the excited harmonics
    if ptree:
        harmonics = array(ptree.get_array_int('harmonics'))
        peak_indices = harmonics * (cycles - ignore_cycles)
    else:
        mx, mn = peakdet(absolute(fft_voltage), mean(absolute(fft_current)))
        peak_indices = int(mx[:, 0])
        mx, mn = peakdet(absolute(fft_voltage), mean(absolute(fft_current)))
        assert peak_indices == mx[:, 0]

    frequency = fft_frequency[peak_indices]
    impedance = fft_voltage[peak_indices] / fft_current[peak_indices]

    return [frequency, impedance]

	def fromEBmodes(cls,fourier_E,fourier_B,angle=3.14*deg):

		"""
		This class method allows to build a shear map specifying its E and B mode components

		:param fourier_E: E mode of the shear map in fourier space
		:type fourier_E: numpy 2D array, must be of type np.complex128 and must have a shape that is appropriate for a real fourier transform, i.e. (N,N/2 + 1); N should be a power of 2

		:param fourier_B: B mode of the shear map in fourier space
		:type fourier_B: numpy 2D array, must be of type np.complex128 and must have a shape that is appropriate for a real fourier transform, i.e. (N,N/2 + 1); N should be a power of 2

		:param angle: Side angle of the real space map in degrees
		:type angle: float.

		:returns: the corresponding ShearMap instance

		:raises: AssertionErrors for inappropriate inputs

		"""

		assert fourier_E.dtype == np.complex128 and fourier_B.dtype == np.complex128
		assert fourier_E.shape[1] == fourier_E.shape[0]/2 + 1
		assert fourier_B.shape[1] == fourier_B.shape[0]/2 + 1
		assert fourier_E.shape == fourier_B.shape

		#Compute frequencies
		lx = rfftfreq(fourier_E.shape[0])
		ly = fftfreq(fourier_E.shape[0])

		#Safety check
		assert len(lx)==fourier_E.shape[1]
		assert len(ly)==fourier_E.shape[0]

		#Compute sines and cosines of rotation angles
		l_squared = lx[np.newaxis,:]**2 + ly[:,np.newaxis]**2
		l_squared[0,0] = 1.0

		sin_2_phi = 2.0 * lx[np.newaxis,:] * ly[:,np.newaxis] / l_squared
		cos_2_phi = (lx[np.newaxis,:]**2 - ly[:,np.newaxis]**2) / l_squared

		sin_2_phi[0,0] = 0.0
		cos_2_phi[0,0] = 0.0

		#Invert E/B modes and find the components of the shear
		ft_data1 = cos_2_phi * fourier_E - sin_2_phi * fourier_B
		ft_data2 = sin_2_phi * fourier_E + cos_2_phi * fourier_B

		#Invert Fourier transforms
		data1 = irfft2(ft_data1)
		data2 = irfft2(ft_data2)

		#Instantiate new shear map class
		new = cls(np.array([data1,data2]),angle)
		setattr(new,"fourier_E",fourier_E)
		setattr(new,"fourier_B",fourier_B)

		return new

def main(input_val):
    currency_pairs = ["BTC_USD_", "LTC_USD_", "LTC_BTC_", "DRK_BTC_",
                      "DOGE_BTC_", "DRK_LTC_", "DOGE_LTC_"]
    
    file_list_uber = []
    
    for pair in currency_pairs:
        file_list_uber.append(pair + "full_ubersampled.csv")
    
    for file in file_list_uber:
        uber_csv = open(file, 'r')
        csv_read = csv.DictReader(uber_csv, delimiter=',',
                                  quoting=csv.QUOTE_NONNUMERIC)
        csv_list = []
        for row in csv_read:
            csv_list.append({'price':float(row['price']),
                             'tmsp':int(row['tmsp'])})            
                    
        first_tmsp = int(csv_list[0]['tmsp'])
        last_tmsp = int(csv_list[-1]['tmsp'])
        delta_tmsp = int(last_tmsp - first_tmsp)
            
        uber_len = len(csv_list)
            
        uber_price_array = np.zeros(uber_len, dtype='float')
            
        for i in range(0,uber_len,1):
            uber_price_array[i] = csv_list[i]['price']
            
        uber_fft = fft.rfft(uber_price_array)
        uber_fft_imag = np.imag(uber_fft)
        uber_fft_real = np.real(uber_fft)
        uber_fft_mag = np.sqrt((uber_fft_real * uber_fft_real) +
                               (uber_fft_imag * uber_fft_imag))
        
        uber_fft_bins = fft.rfftfreq(uber_len, d=100)
        
        uber_csv.close()
                
        field_names = ['freq_bin', 'magnitude']
        fft_csv_name = file[0:-4] + "_fft.csv"
        
        fft_csv = open(fft_csv_name, 'w', newline = '')
                
        csvwriter = csv.DictWriter(fft_csv, delimiter=',',fieldnames=field_names,
                                   quoting=csv.QUOTE_NONNUMERIC)
        csvwriter.writeheader()
        #csvwriter.writerows(full_conditioned)
                
        for i in range(0,len(uber_fft_real),1):
            dict_d = {'freq_bin':uber_fft_bins[i],
                      'magnitude':uber_fft_mag[i]}
            csvwriter.writerow(dict_d)
                    
        fft_csv.close()        

def freq_analyze(data, noise, dt):
    plt.figure()
    plt.plot(data)
    plt.title('data to analyze')
    
    plt.figure()
    plt.title('noise to analyze')
    plt.plot(noise)
    
    data_FFT=rfft(data)
    data_freqs=rfftfreq(len(data), dt)
    noise_FFT=rfft(noise)
    noise_freqs=rfftfreq(len(noise), dt)
    
    plt.figure()
    plt.plot(data_freqs, data_FFT, label='data FFT')
    plt.plot(noise_freqs, noise_FFT, label='noise FFT')
    plt.title('FFT of noise and data')
    plt.yscale('log')
    plt.show()

def getSTFT(signal, sample_spacing=1, slidingWindowSize=20, stepSize=10) :
    """
    Get short-time fourier transform with slidingWindowSize=20 (number of sample in window) and stepSize=10 (number of samle to sample from window to another)
    """
    segmentedSignal=segmentSignal(signal,slidingWindowSize,stepSize)
    segmentedSignalSTFT=[]
    for segment in segmentedSignal : segmentedSignalSTFT.append(rfft(segment))
    spectrogram=np.array(segmentedSignalSTFT)
    t=np.array([stepSize*i for i in xrange(len(segmentedSignal))])
    f=np.array(rfftfreq(slidingWindowSize,d=sample_spacing))
    return t,f,spectrogram

 def __init__(self, args):
     self.args = args
     # Determine the frequencies the dfft calculates at.
     self.freq = fft.rfftfreq(n=args.fft_size,
                              d=(1 / args.sample_rate))
     # Setup the display:
     self.fig = plt.figure()
     self.ax = plt.axes(xlim=args.xlims, ylim=args.ylims)
     self.ax.set_xscale('log', basex=10)
     self.line, = self.ax.plot([], [])
     plt.ion()

    def __init__(self, **kwargs):
        super(SpectralFrameWidget, self).__init__()

        # default parameters if not specified
        for name, val in self.defaults.items():
            setattr(self, name, kwargs.get(name, val))

        self.setLabel('left', 'Amplitude', units='dB')
        self.setLabel('bottom', text='Frequency', units='kHz')
        self.getPlotItem().setRange(yRange=[-50, 40])
        self.spectrum_curve = self.plot([])
        self.freqs = rfftfreq(self.nfft) * self.fs / 1000.0
        self.frame_size = get_frame_size(self.fs, self.frame_size_ms)
        self.window = np.hamming(self.frame_size)

def simpleSpecs(t, sig, amponly=False):
    '''
    Creates simple amplitude and phase spectrums of real functions using numpy.fft; phase is in degrees
    :param t: time (1D numpy array)
    :param sig: signal (1D numpy array)
    :param amponly: outputs only the amplitude spectrum (bool)
    :return: frequency (1D numpy array), A(f), phi(f)
    '''
    n = sig.size
    fouSig = fft.rfft(sig) / n # DFT of the signal normalized
    sample_int = t[n-1] - t[n-2] # sample interval
    fouFrq = fft.rfftfreq(n, d=sample_int)
    fouAmp = sqrt(real(fouSig)**2 + imag(fouSig)**2)
    fouPhi = angle(fouSig, deg=True)
    if amponly is True:
        return fouFrq, fouAmp
    else:
        return fouFrq, fouAmp, fouPhi

	def init_plot(self):
		x = fft.rfftfreq(self.CHUNK_SIZE)*self.SAMPLE_RATE
		x_max = x[-1]
		self.init_y = linspace(0, x_max, len(x))
		y = self.init_y
		source = ColumnDataSource(data=dict(x=x, y=y))
		
		# TODO: range and size (toolbar), maybe could be user settings
		plot = Figure(plot_height=400, plot_width=800, title="freq anal",
					  tools="crosshair,pan,reset,resize,save,wheel_zoom",
					  x_range=[0, x_max], y_range=[0, 15])
					  
		rad = x_max/float(len(x))
		data = plot.circle('x', 'y', source=source, line_color=None, radius=rad)
		self.data_source = data.data_source
		
		# TODO: maybe not the best place
		curdoc().add_root(plot)

def testResample():
  s, rate = ReadAndConvert(os.path.join('testdata', 'foo_s80_p95.wav'))
  s = scipy.signal.resample(s, 4*len(s))
  rate = 4*rate
  ConvertAndWrite(s, rate, os.path.join('testout', 'foo_resampled.wav'))
  s = s[0 : 15*2048]
  nfft, nskip = 1024, 512
  Sws = STFT(s, nfft, nskip)
  nrfft, nrskip = 4096, 2048
  target_pitch = 880. #440.0 * 2 ** (7.0 / 12)

  Rws = []
  for nw, Sw in enumerate(Sws):
    freqs = rfftfreq(nfft, 1.0/rate)
    epitch = EstimatePitch(Sw, rate)

    if epitch == 0.0: 
      p = 1
    else:
      p = target_pitch / epitch
      print '%d: epitch=%f p=%f' % (nw, epitch, p)

    PlotPitches(np.abs(Sw), rate, name='testResampleFrame%dInMag' % nw,
                title='Spectrum of input frame %d (%f)' % (nw, epitch))
    PlotPitches(np.angle(Sw), rate, name='testResampleFrame%dInPhase' % nw,
                title='Phase of input frame %d' % nw)

    Rw = Resample(Sw, nfft, nrfft, p, rate)

    Rws.append(Rw)

    erpitch = EstimatePitch(Rw, rate)
    PlotPitches(np.abs(Rw), rate, name='testResampleFrame%dOutMag' % nw,
                title='Spectrum of output frame %d (%f)' % (nw, erpitch))
    PlotPitches(np.angle(Rw), rate, name='testResampleFrame%dOutPhase' % nw,
                title='Phase of output frame %d' % nw)
  r = STIFT(Rws, nrfft, nrskip)
  r = scipy.signal.resample(r, int(np.floor(r.size/2)))
  rate = rate/2

  ConvertAndWrite(
      r, rate,
      os.path.join('testout', 'testResample.wav'))

def real_to_fourier(f, x, y, z):
    # Check that real-space grid spacing is all equal
    if not (_is_evenly_spaced(x) and _is_evenly_spaced(y) and _is_evenly_spaced(z)):
        raise ValueError('Sample points in real space are not evenly spaced.')
    dx = x[1]-x[0]  # Grid spacing
    dy = y[1]-y[0]
    dz = z[1]-z[0]

    ftrans  = ft.rfftn(f)

    # Wavenumber arrays
    kx        = 2*np.pi * ft.fftfreq(x.size, d=dx)
    ky        = 2*np.pi * ft.fftfreq(y.size, d=dy)
    kz        = 2*np.pi * ft.rfftfreq(z.size, d=dz) # Only last axis is halved in length when using numpy.fft.rfftn()

    # Normalize (convert DFT to continuous FT)
    ftrans *= dx*dy*dz

    return ftrans, kx, ky, kz

    def PlotPSD(self):
        colors = ['r']
        dt = self.Dt
        Ry = self.measurements[0][0]
        mean = np.average(Ry)
        val = np.array( [v - mean for v in Ry] )
        freq = fft.rfftfreq(len(val), dt)
        tr = fft.rfft(val)
        tr_vect = [np.real(i) for i in tr]

        psd = [(np.real(i)**2 + np.imag(i)**2)**0.5 for i in tr]

        plt.clf()
        plt.figure("psd")
        plt.plot(freq[1:],tr_vect[1:],colors[0])
        plt.xlim(0, 10.0)
        plt.xlabel(r"$Frequency \, (Hz)$")
        plt.ylabel(r"$PSD$")
        plt.savefig("psd.eps")

    def fft(self, nfft=None):
        """Compute the one-dimensional discrete Fourier transform of
        this `TimeSeries`.

        Parameters
        ----------
        nfft : `int`, optional
            length of the desired Fourier transform.
            Input will be cropped or padded to match the desired length.
            If nfft is not given, the length of the `TimeSeries`
            will be used

        Returns
        -------
        out : :class:`~gwpy.spectrum.Spectrum`
            the normalised, complex-valued FFT `Spectrum`.

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

        Notes
        -----
        This method, in constrast to the :meth:`numpy.fft.rfft` method
        it calls, applies the necessary normalisation such that the
        amplitude of the output :class:`~gwpy.spectrum.Spectrum` is
        correct.
        """
        from ..spectrum import Spectrum
        if nfft is None:
            nfft = self.size
        dft = npfft.rfft(self.value, n=nfft) / nfft
        dft[1:] *= 2.0
        new = Spectrum(dft, epoch=self.epoch, channel=self.channel,
                       unit=self.unit)
        new.frequencies = npfft.rfftfreq(self.size, d=self.dx.value)
        return new