def phi(entry):
    # Likelihood fonction
    A = MesFonctions.matrice_A(entry[0:size].reshape((NbLigne, NbColonne)), methode="3")
    x_ = (
        floue.reshape(size) - MesFonctions.conv_matrices(A, entry[size:2 * size].reshape((NbLigne, NbColonne)),
                                                         "vecteur",
                                                         "coin"))
    likelihood = np.dot(np.dot(x_, C_noise_inv), x_.transpose())
    # Regularization term
    image_ = entry[0:size]  # /sum(sum(im))
    PSF_ = entry[size:2 * size]  # /sum(sum(PSF))
    if regularization == "spectral":
        fy_ = np.fft(image_)
        y_ = 0
        for i in range(0, size):
            y_ = y_ + i * i * abs(fy_[i]) * abs(fy_[i])
        regularization_im = y_
        fz_ = np.fft(PSF_)
        z_ = 0
        for i in range(0, size):
            z_ = z_ + i * i * abs(fz_[i]) * abs(fz_[i])
        regularization_PSF = z_
    elif regularization == "Tikhonov":
        y_ = sum(image_ * image_)
        regularization_im = math.sqrt(y_)
        z_ = sum(PSF_ * PSF_)
        regularization_PSF = math.sqrt(z_)

    t_ = likelihood + regularization_im + regularization_PSF
    return t_

	def process2ch(self,inputbuffers,timestamp):
		if not self.update() : return None

		fftsize = self.m_blockSize
		
		audioSamples0 = inputbuffers[0]
		audioSamples1 = inputbuffers[1]

		complexSpectrum0 = fft(self.window*audioSamples0,fftsize)
		complexSpectrum1 = fft(self.window*audioSamples1,fftsize)
		
		magnitudeSpectrum0 = abs(complexSpectrum0)[0:fftsize/2] / (fftsize/2)
		magnitudeSpectrum1 = abs(complexSpectrum1)[0:fftsize/2] / (fftsize/2)
		
		# do the computation
		melSpectrum0 = self.warpSpectrum(magnitudeSpectrum0)
		melCepstrum0 = self.getMFCCs(melSpectrum0,cn=True)
		melSpectrum1 = self.warpSpectrum(magnitudeSpectrum1)
		melCepstrum1 = self.getMFCCs(melSpectrum1,cn=True)
		
		outputs = FeatureSet()
		outputs[0] = Feature(hstack((melCepstrum1[self.cnull:],melCepstrum0[self.cnull:])))
		outputs[1] = Feature(hstack((melSpectrum1,melSpectrum0)))
		
		return outputs

    def calculate_FROG(self):
        
        self.FROGxmin   = self.t[0]
        self.FROGxmax   = self.t[-1]
        self.FROGdeltax = (self.t[-1]-self.t[0])/self.NT
        self.FROGymin   = self.frequencies[0]
        self.FROGymax   = self.frequencies[-1]
        self.FROGdeltay = (self.frequencies[-1]-self.frequencies[0])/self.NT
        return
#negative delay
        for i in xrange(self.NT/2):
            self.field     =  zeros((self.NT))
            self.field[i:] = (self.ElectricField[:self.NT-i]*self.ElectricField[i:])**2
            self.field_fft = fft(self.field)       
            self.field_fft[:self.NT/2] = zeros((self.NT/2))   #no negative freq
            self.FROG[:,self.NT/2-i] = abs(self.field_fft[self.NT/2:])**2#fftshift(abs(self.field_fft[self.NT/2:])**2)
            
#positive delay
        for i in xrange(self.NT/2):
            self.field     = zeros((self.NT))
            self.field[i:] = (self.ElectricField[i:]*self.ElectricField[:self.NT-i])**2
            self.field_fft = fft(self.field)  
            self.field_fft[:self.NT/2] = zeros((self.NT/2))      #no negative freq
            self.FROG[:,self.NT/2+i] = abs(self.field_fft[self.NT/2:])**2#fftshift(abs(self.field_fft[self.NT/2:])**2)
            
        self.FROG /= ma.max(abs(self.FROG))

    def calculate_FROG(self):
        
        self.FROGxmin   = self.t[0]
        self.FROGxmax   = self.t[-1]
        self.FROGymin   = self.get_frequencies()[0]
        self.FROGymax   = self.get_frequencies()[-1]
        
        ElectricField = real(self.ElectricField)
        # complex field
        
#negative delay
        for i in xrange(self.NT/2):
            self.field     =  zeros((self.NT))
            self.field[i:] = (ElectricField[:self.NT-i]*ElectricField[i:])**2
            self.field_fft = fft(self.field)       
            self.FROG[:,self.NT/2-i] = fftshift(abs(self.field_fft)**2)
            
#positive delay
        for i in xrange(self.NT/2):
            self.field     = zeros((self.NT))
            self.field[i:] = (ElectricField[i:]*ElectricField[:self.NT-i])**2
            self.field_fft = fft(self.field)  
            self.FROG[:,self.NT/2+i] = fftshift(abs(self.field_fft)**2)
            
        self.FROG /= ma.max(abs(self.FROG))

    def myFFT(self,x):
        tt,xx,uu = self.Traj(x,U=1,fast=0)
        n,gn=len(tt),len(uu)
        grsU = [TGraph(n,tt,uu[i]) for i in range(5)]         
        ffts = [fft(uu[i]).real for i in range(5)]
        ffts2= [fft(uu[i]).imag for i in range(5)]
        fftUt = [array('d') for i in range(5)]
        fftUt2 = [array('d') for i in range(5)]
        for i in range(5):
            for j in range(n):
                fftUt[i].append(ffts[i][j]**2+ffts2[i][j]**2)
                fftUt2[i].append(ffts[i][j])
        fftUr = [TGraph(n,tt,fftUt[i]) for i in range(5)]         
        fftUi = [TGraph(n,tt,fftUt2[i]) for i in range(5)]         
        c1 = TCanvas()
        c1.Divide(3,2)
        for i in range(5):
            c1.cd(i+1)
            fftUi[i].Draw("APE")
 
        fftsx1 = [fft(xx[i]).real for i in range(3)]
        fftsx2= [fft(xx[i]).imag for i in range(3)]
        fftsxA = [array('d') for i in range(3)]
        for i in range(3):
            for j in range(n):
                fftsxA[i].append(fftsx1[i][j])
        fftxr = [TGraph(n,tt,fftsxA[i]) for i in range(3)]         
        
        c2 = TCanvas()
        c2.Divide(3)
        
        for i in range(3):
            c2.cd(i+1)
            fftxr[i].Draw("APE")
        zwom = input("numbers continue")

def length_shift(signalin, tscale):
  L = len(signalin)
  # signal blocks for processing and output
  phi  = zeros(N)
  out = zeros(N, dtype=complex)
  sigout = zeros(L/tscale+N)

  # max input amp, window
  amp = max(signalin)
  win = hanning(N)
  p = 0
  pp = 0
  while p < L-(N+H):
    print p, ',', (L-N-H)

    # take the spectra of two consecutive windows
    p1 = int(p)
    spec1 =  fft(win*signalin[p1:p1+N])
    spec2 =  fft(win*signalin[p1+H:p1+N+H])
    # take their phase difference and integrate
    phi += (angle(spec2) - angle(spec1))
    # bring the phase back to between pi and -pi
    for i in phi:
      while i < -pi: i += 2*pi
      while i >= pi: i -= 2*pi
    out.real, out.imag = cos(phi), sin(phi)
    # inverse FFT and overlap-add
    sigout[pp:pp+N] += (win*ifft(abs(spec2)*out)).real
    pp += H
    p += H*tscale
  return array(amp * sigout / max(sigout))

def fft_sheet(f, t):
    # Transform along y
    fk = fft(f, axis=0)
    # Phase shift (to make periodic in x)
    fk *= exp(1j * S * t * ky * x)
    # Transform along x
    return fft(fk, axis=1)

def fft_covariance(sample):
    '''
        returns the covariance matrix and index of each item for a dimension-list of sample-lists
        [d0, d1, ...,dd] where d0 is the first dimension containing n samples.
        n.b. numpy 1.7rc1 has numpy.pad to pad array.  currently running 1.6.2 (released version)
        so imported rc1 version.
    '''
        
    A = numpy.array(sample)
    d, n = A.shape
    Q = numpy.ones(1) if(d==2) else numpy.ones((d, d))
    idx = [[1]] if(d==2) else [list([1.] * d) for i in range(d)]
    nm1 = float(n - 1)
    #prep for fft
    A =  (A.transpose() - numpy.mean(A, axis=1)).transpose()
    A = npp.pad(A,(0,len(pad(n))), 'constant', constant_values=(0,0))[:d, : ]
    #2x2return
    if (d == 2):
        corr = ifft(fft(A[0, :]) * np.conjugate(fft(A[1, :])))
        maxcorr = getmax(corr)
        Q[0] =  maxcorr[0] / nm1
        idx[0] = maxcorr[1]
    else:
        for i in range(d):
            for j in range(i, d):
                corr = ifft(fft(A[i, :]) * np.conjugate(fft(A[j, : ])))
                maxcorr = getmax(corr)
                Q[i][j] = maxcorr[0] / nm1 
                idx[i][j] = maxcorr[1]
                if i != j:
                    Q[j][i] = Q[i][j]
    return Q, idx

def FFT_Correlation(x,y):
    """
    FFT-based correlation, much faster than numpy autocorr.
    x and y are row-based vectors of arbitrary lengths.
    This is a vectorized implementation of O(N*log(N)) flops.
    """

    lengthx = x.shape[0]
    lengthy = y.shape[0]

    x = np.reshape(x,(1,lengthx))
    y = np.reshape(y,(1,lengthy))

    length = np.array([lengthx, lengthy]).min()
    
    x = x[:length]
    y = y[:length]
    
    fftx = fft(x, 2 * length - 1, axis=1) #pad with zeros
    ffty = fft(y, 2 * length - 1, axis=1)

    corr_xy = fft.ifft(fftx * np.conjugate(ffty), axis=1)
    corr_xy = np.real(fft.fftshift(corr_xy, axes=1)) #should be no imaginary part

    corr_yx = fft.ifft(ffty * np.conjugate(fftx), axis=1)
    corr_yx = np.real(fft.fftshift(corr_yx, axes=1))

    corr = 0.5 * (corr_xy[:,length:] + corr_yx[:,length:]) / range(1,length)[::-1]
    return np.reshape(corr,corr.shape[1])

 def diff_spectrum(self, ref_mic):
     ref_signal = self.microphones[ref_mic[0]][ref_mic[1]].signal
     min_sig_len = min(ref_signal.shape[0], self.signal.shape[0])
     wnd = 0.54 - 0.46*numpy.cos(2*pi*numpy.arange(min_sig_len)/min_sig_len)
     fft_sum  = fftshift(fft(self.signal[0:min_sig_len]*wnd))/min_sig_len
     fft_ref = fftshift(fft(ref_signal[0:min_sig_len]*wnd))/min_sig_len
     fft_rel  = (fft_sum / fft_ref).real
     fft_freq = numpy.linspace(-self.samplerate/2, self.samplerate/2, min_sig_len);
     return {"x":fft_freq, "y":fft_rel}

def fftconv(x, y):
    """ Convolution of x and y using the FFT convolution theorem. """
    N = len(x)
    n = int(2 ** np.ceil(np.log2(N))) + 1
    X, Y, x_y = fft(x, n), fft(y, n), []
    for i in range(n):
        x_y.append(X[i] * Y[i])

    # Returns the inverse Fourier transform with padding correction
    return fft.ifft(x_y)[4:N+4]

def conv(x,y):
    x1=np.zeros(x.size) 
    x1[0:x.size]=x  #added zeros

    y1=np.zeros(y.size)
    y1[0:y.size]=y #added zeros
    x1ft=np.fft(x1)
    y1ft=np.fft(y1)
    vect1=np.real(np.ifft(x1ft*y1ft))
    return vect1[0:x.size]

def my_fft(sample, inverse=False):
    n = len(sample)
    if n == 1:
        return sample
    else:
        Feven = fft([sample[i] for i in xrange(0, n, 2)])
        Fodd = fft([sample[i] for i in xrange(1, n, 2)])

        combined = [0] * n
        for m in xrange(n/2):
            combined[m] = Feven[m] + omega(n, -m, inverse) * Fodd[m]
            combined[m + n/2] = Feven[m] - omega(n, -m, inverse) * Fodd[m]

        return combined

 def solve(self, x, rho):
     Nx, = rho.shape
     k = fftfreq(Nx) * 2*np.pi*Nx
     k[0] = 1.0
     rhohat = fft(rho)
     rhobar = rhohat[0] / rho.size
     phihat = rhohat / -k**2 + fft(0.5 * rhobar * x**2)
     gphhat = 1.j * k * phihat
     phi = ifft(phihat).real
     gph = ifft(1.j * k * phihat).real
     soln = np.zeros(rho.shape + (4,))
     soln[:,0] = phi
     soln[:,1] = gph
     return soln

def showResponse(x, f, a, w0, Q):
	fmax = 200.0
	
	# smooth it out with averaging:
	df = 1.0
	fatf = arange(0,fmax,df)
	n = len(fatf)
	fatX = zeros(n)

	print '\nfft\'ing to find the frequency response...'
	X = fft(x)
	for i in range(n):
		fslice = where(abs(f-fatf[i])<df/2)[0]
		fatX[i] = mean(abs(X[fslice]))

	norm = mean(fatX[0:10])
	fatX = fatX/norm

	# plot the response with the theoretical curve overlaid:
	plot(fatf, fatX, 'k.', label = 'simulated data')
	title('response of damping system')
	xlabel('frequency (Hz)')
	ylabel('abs(fft(x)) normalized to low f')

	Omega = 2*pi*fatf/w0
	zeta = 1.0/2.0/Q
	theory = 1/sqrt( (1-Omega**2)**2 + (2*zeta*Omega)**2 )
	
	semilogy(fatf, theory, label = 'theory')
	legend()
	show()

def find_frequency(x,y):		# Returns the fundamental frequency using FFT
	tx,ty = fft(y, x[1]-x[0])
	index = find_peak(ty)
	if index == 0:
		return None
	else:
		return tx[index]

def fft_squid(path):
    squid_model = model.NoisySquid(path)
    squid_model.getContext().setClock(3, 20e-3, 0)
    squid_model.currentInjection.min = 0.0e-6
    squid_model.currentInjection.max = 0.1e-6
    squid_model.getContext().reset()
    squid_model.run(200e-3)
    signal = array(squid_model.vmTable)
    squid_model.save_all_plots()
#    savetxt('vm.plot', signal)
    transform = fft(signal)
    xx = array(range(len(transform)))
    ax1 = subplot(121)
    #plot(xx / 10.0 , signal * 1e3, 'r')
    savetxt("t_series.plot", xx/10.0)
#     ax1 = axis([0, 0, 0.5, 0.8])
#     xlabel("time (ms)")
#     ylabel("membrane potential (mV)")
#     ax2 = twinx()
#     plot(xx / 10.0 , array(squid_model.iInjectTable) * 1e9, 'b--')
#     ylabel("injection current (nA)")
#     ax2.yaxis.tick_right()
#     title("(A)", fontname="Times New Roman", fontsize=10, fontweight="bold")
#     show()
#     savefig("Vm.jpg", dpi=600)
    
    savetxt('fftreal.plot', transform.real, fmt="%13.12G")
    savetxt('fftimag.plot', transform.imag, fmt="%13.12G")

 def fourier_projection(self, density):
     
     ft_density = np.fft(density)
     
     # -m wraps around array, but that should be OK....
     # --> self._A_ell_expt[l] is shape (2l+1) x n_q  = (m x q)
     #     ft_coefficients[:,l,:] is shape n_q x (2l+1) = (q x m) (tranposed later)
     ft_coefficients = self._sph_harm_projector(ft_density)
     
     
     # zero out the array, use it to store the next iter
     ft_density[:,:,:] = 0.0 + 1j * 0.0
     
     for l in range(self.order_cutoff):
         
         A_ell_model = ft_coefficients[:,l,:].T # (2l+1) x n_q  = (m x q)
     
         # find U that rotates the experimental vectors into as close
         # agreement as possible as the model vectors
         U = math2.kabsch(self._A_ell_expt[l], A_ell_model)
 
         # update: k --> k+1
         A_ell_prime = np.dot(self._A_ell_expt[l], U)
         ft_density_prime += self._sph_harm_projector.expand_sph_harm_order(A_ell_prime, l)
         
     updated_density = np.ifft(ft_density)
         
     return updated_density

    def processBuffer(self, bounds):
        self.param1 = bounds.height/65536.0
        self.param2 = bounds.height/2.0

        if(self.stop==False):

            Fs = 48000
            nfft= 65536
            self.newest_buffer=self.newest_buffer[0:256]
            self.fftx = fft(self.newest_buffer, 256,-1)

            self.fftx=self.fftx[0:self.freq_range*2]
            self.draw_interval=bounds.width/(self.freq_range*2.)

            NumUniquePts = ceil((nfft+1)/2)
            self.buffers=abs(self.fftx)*0.02
            self.y_mag_bias_multiplier=0.1
            self.scaleX = "hz"
            self.scaleY = ""

        if(len(self.buffers)==0):
            return False

        # Scaling the values
        val = []
        for i in self.buffers:
            temp_val_float = float(self.param1*i*self.y_mag) + self.y_mag_bias_multiplier * self.param2

            if(temp_val_float >= bounds.height):
                temp_val_float = bounds.height-25
            if(temp_val_float <= 0):
                temp_val_float = 25
            val.append( temp_val_float )

        self.peaks = val

def dst(x,axis=-1):
    """Discrete Sine Transform (DST-I)

    Implemented using 2(N+1)-point FFT
    xsym = r_[0,x,0,-x[::-1]]
    DST = (-imag(fft(xsym))/2)[1:(N+1)]

    adjusted to work over an arbitrary axis for entire n-dim array
    """
    n = len(x.shape)
    N = x.shape[axis]
    slices = [None]*3
    for k in range(3):
        slices[k] = []
        for j in range(n):
            slices[k].append(slice(None))
    newshape = list(x.shape)
    newshape[axis] = 2*(N+1)
    xtilde = np.zeros(newshape,np.float64)
    slices[0][axis] = slice(1,N+1)
    slices[1][axis] = slice(N+2,None)
    slices[2][axis] = slice(None,None,-1)
    for k in range(3):
        slices[k] = tuple(slices[k])
    xtilde[slices[0]] = x
    xtilde[slices[1]] = -x[slices[2]]
    Xt = np.fft(xtilde,axis=axis)
    return (-np.imag(Xt)/2)[slices[0]]