def MoyalPropagation(W):
    """
    Propagate wigner function W by the Moyal equation.
    This function is used to verify that the obtained wigner functions
    are steady state solutions of the Moyal equation.
    """
    # Make a copy
    W = np.copy(W)

    dt = 0.005 # time increment
    TIterSteps = 2000

    # Pre-calculate exp
    expIV = np.exp(-1j*dt*(V(X - 0.5*Theta) - V(X + 0.5*Theta)))
    expIK = np.exp(-1j*dt*(K(P + 0.5*Lambda) - K(P - 0.5*Lambda)))

    for _ in xrange(TIterSteps):

        # p x -> theta x
        W = fftpack.fft(W, axis=0, overwrite_x=True)
        W *= expIV
        # theta x  ->  p x
        W = fftpack.ifft(W, axis=0, overwrite_x=True)

        # p x  ->  p lambda
        W = fftpack.fft(W, axis=1, overwrite_x=True)
        W *= expIK
        # p lambda  ->  p x
        W = fftpack.ifft(W, axis=1, overwrite_x=True)

        # normalization
        W /= W.real.sum()*dX*dP

    return fftpack.fftshift(W.real)

def highpass_filter(data, width):
    """Highpass filter on *width* scales using blackman window.

    Finite impulse response filter *that discards invalid data* at the ends.

    """

    ntime = data.shape[-1]

    # Blackman FWHM factor.
    window_width = int(width / 0.4054785)

    if window_width % 2:
        window_width += 1

    window = np.zeros(ntime, dtype=np.float32)
    window_core = signal.blackman(window_width, sym=True)
    window_core = -window_core / np.sum(window_core)
    window_core[window_width // 2] += 1
    window[:window_width] = window_core
    window_fft = fftpack.fft(window)

    ntime_out = data.shape[-1] - window_width + 1
    out_shape = data.shape[:-1] + (ntime_out,)
    out = np.empty(out_shape, data.dtype)

    for ii in range(data.shape[0]):
        d_fft = fftpack.fft(data[ii])
        d_fft *= window_fft
        d_lpf = fftpack.ifft(d_fft)
        out[ii] = d_lpf[-ntime_out:].real

    return out

def convolve_scalogram(ana, wf, sampling_rate,optimize_fft):
    n = wf.shape[0]
    sig  = ana.magnitude
    ana_sr=ana.sampling_rate.rescale('Hz').magnitude
    if optimize_fft:
        sig=sig-sig.mean() # Remove mean before padding
        nfft=int(2**np.ceil(np.log(sig.size)/np.log(2)))
        sig=np.r_[sig,np.zeros(nfft-sig.size)] # pad signal with 0 to a power of 2 length
        sig=resample(sig,int(sig.size*sampling_rate/ana_sr)) # resample in time domain 
        sigf=fftpack.fft(sig,n) # Compute fft with a power of 2 length        
    else:        
        sigf=fftpack.fft(sig)
        # subsampling in fft domain (attention factor)
        factor = (sampling_rate/ana.sampling_rate).simplified.magnitude
        x=(n-1)//2
        if np.mod(n,2)==0:
            sigf = np.concatenate([sigf[0:x+2],  sigf[-x:]])*factor
        else:
            sigf = np.concatenate([sigf[0:x+1],  sigf[-x:]])*factor
            
    # windowing ???
    #win = fftpack.ifftshift(np.hamming(n))
    #sigf *= win
    
    # Convolve (mult. in Fourier space)
    wt_tmp=fftpack.ifft(sigf[:,np.newaxis]*wf,axis=0)
    # and shift
    wt = fftpack.fftshift(wt_tmp,axes=[0])
    return wt

def fcglt(A): # Modal Coefficients to Lobatto Nodal
    """
    Fast Chebyshev-Gauss-Lobatto transformation from 
    Chebyshev expansion coefficients (modal) to point 
    space values (nodal). If I=numpy.identity(n), then
    T=chebyshev.fcglt(I) will be the Chebyshev 
    Vandermonde matrix on the Lobatto nodes
    """
    size = A.shape
    m = size[0]
    k = m-2-np.arange(m-2)

    if len(size) == 2: # Multiple vectors
        V = np.vstack((2*A[0,:],A[1:m-1,:],2*A[m-1,:],A[k,:]))
        F = fft(V, n=None, axis=0)
        B = 0.5*F[0:m,:]
    else:  # Single vector
        V = np.hstack((2*A[0],A[1:m-1],2*A[m-1],A[k]))
        F = fft(V, n=None)
        B = 0.5*F[0:m]

    if A.dtype!='complex':
        return np.real(B)
    else:
        return B

    def run(self):
        if self.filter_on == False:
            f = lambda x:random.random()+self.amp*np.sin(x)
            x = np.linspace(0, 10)

            fft_feature = fft.fft(map(f, x))

            ifft_feature = fft.ifft(fft_feature)

            self.newSample.emit(map(f, x))
            self.newSamplefft.emit(list(abs(fft_feature)))
            self.newSampleifft.emit(list(ifft_feature))
        elif self.filter_on == True:
            f = lambda x:random.random()+self.amp*np.sin(x)
            x = np.linspace(0, 10)

            fft_feature = fft.fft(map(f, x))
            mean = np.average(abs(fft_feature))
            fft_feature_filter = fft_feature
            for i in range(len(fft_feature)):
                if abs(fft_feature[i]) >= mean:
                    fft_feature_filter[i] = abs(fft_feature[i])
                else:
                    fft_feature_filter[i] = 0

            ifft_feature = fft.ifft(fft_feature_filter)

            self.newSample.emit(map(f, x))
            self.newSamplefft.emit(list(fft_feature_filter))
            self.newSampleifft.emit(list(ifft_feature))

            self.filter_on = False
        else:
            pass

def make_audio_analysis_plots(infile, prefix='temp', make_plots=True,
                              do_fft=True, fft_sum=None):
    ''' create frequency plot '''
    import numpy as np
    from scipy import fftpack
    from scipy.io import wavfile
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as pl

    if not os.path.exists(infile):
        return -1

    try:
        rate, data = wavfile.read(infile)
    except ValueError:
        print('error reading wav file')
        return -1
    dt_ = 1./rate
    time_ = dt_ * data.shape[0]
    tvec = np.arange(0, time_, dt_)
    sig0 = data[:, 0]
    sig1 = data[:, 1]
    if not tvec.shape == sig0.shape == sig1.shape:
        return -1
    if not do_fft:
        fft_sum_ = float(np.sum(np.abs(sig0)))
        if hasattr(fft_sum, 'value'):
            fft_sum.value = fft_sum_
        return fft_sum_
    if make_plots:
        pl.clf()
        pl.plot(tvec, sig0)
        pl.plot(tvec, sig1)
        xtickarray = range(0, 12, 2)
        pl.xticks(xtickarray, ['%d s' % x for x in xtickarray])
        pl.savefig('%s/%s_time.png' % (HOMEDIR, prefix))
        pl.clf()
    samp_freq0 = fftpack.fftfreq(sig0.size, d=dt_)
    sig_fft0 = fftpack.fft(sig0)
    samp_freq1 = fftpack.fftfreq(sig1.size, d=dt_)
    sig_fft1 = fftpack.fft(sig1)
    if make_plots:
        pl.clf()
        pl.plot(np.log(np.abs(samp_freq0)+1e-9), np.abs(sig_fft0))
        pl.plot(np.log(np.abs(samp_freq1)+1e-9), np.abs(sig_fft1))
        pl.xlim(np.log(10), np.log(40e3))
        xtickarray = np.log(np.array([20, 1e2, 3e2, 1e3, 3e3, 10e3, 30e3]))
        pl.xticks(xtickarray, ['20Hz', '100Hz', '300Hz', '1kHz', '3kHz',
                               '10kHz', '30kHz'])
        pl.savefig('%s/%s_fft.png' % (HOMEDIR, prefix))
        pl.clf()

        run_command('mv %s/%s_time.png %s/%s_fft.png %s/public_html/videos/'
                    % (HOMEDIR, prefix, HOMEDIR, prefix, HOMEDIR))

    fft_sum_ = float(np.sum(np.abs(sig_fft0)))
    if hasattr(fft_sum, 'value'):
        fft_sum.value = fft_sum_
    return fft_sum_

    def fitTrace(self,kwidth=10,porder=3,cwidth=30,pad=False):
        sh = self.sh
        xr1 = (0,sh[1])
        xrs = [xr1]

        polys = []
        for xr in xrs:
            xindex = np.arange(xr[0],xr[1])
            kernel = np.median(self.image[int(sh[0]/2-kwidth):int(sh[0]/2+kwidth),xindex],0)
                
            centroids = []
            totals = []
            for i in np.arange(sh[0]):
                row = self.image[i,xindex]
                row_med = np.median(row)
                    
                total = np.abs((row-row_med).sum())
                cc = fp.ifft(fp.fft(kernel)*np.conj(fp.fft(row-row_med)))
                cc_sh = fp.fftshift(cc)
                centroid = helpers.calc_centroid(cc_sh,cwidth=cwidth).real - xindex.shape[0]/2.
                centroids.append(centroid)
                totals.append(total)

            centroids = np.array(centroids)
        
            yindex = np.arange(sh[0])
            gsubs = np.where((np.isnan(centroids)==False))

            centroids[gsubs] = median_filter(centroids[gsubs],size=20)
            coeffs = np.polyfit(yindex[gsubs],centroids[gsubs],porder)

            poly = np.poly1d(coeffs)
            polys.append(poly)
        return xrs,polys

def freq_filter(f, dt, cutoff_freq, convention='math'):
    """
    A digital filter that filters frequency below a cutoff frequency

    Parameters
    ----------
    f : time signal
    dt : sampling period
    nu_cutoff : cutoff frequency

    Returns
    -------
    The filtered time signal
    """

    if convention == 'math':
        f_freq = fft(f)
    elif convention == 'physics':
        f_freq = ifft(f)
    Ns = np.size(f)
    freqs = fftfreq(Ns, dt)

    # filtering operation
    f_freq[np.where(np.abs(freqs) > cutoff_freq)] = 0
    # go back to time domain
    if convention == 'math':
        f_filter_time = ifft(f_freq)
    elif convention == 'physics':
        f_filter_time = fft(f_freq)

    return f_filter_time

def calc_QW(n, k, kk, kw,  Q, W, useFFT):
    """
    Convolution coefficient

    Args:
        n: Order of the coefficients

        k: Index of the bus

        C: Voltage coefficients (Ncoeff x nbus elements)

    Output:
        Convolution coefficient of order n for the bus k
    """

    if useFFT:
        a = fftpack.fft(Q[:, kk])
        b = fftpack.fft(conj(W[:, kw]))
        e = fftpack.ifft(a * b)
        result = e[n]
    else:
        result = complex_type(0)
        for l in range(n):
            result += Q[l, kk] * conj(W[n-l, kw])

    return result

def l0_gradient_minimization_1d(I, lmd, beta_max, beta_rate=2.0, max_iter=30, return_history=False):
    S = np.array(I).ravel()

    # prepare FFT
    F_I = fft(S)
    F_denom = np.abs(psf2otf([-1, 1], S.shape[0]))**2.0

    # optimization
    S_history = [S]
    beta = lmd*2.0
    hp = np.zeros_like(S)
    for i in range(max_iter):
        # with S, solve for hp in Eq. (12)
        hp = circulant_dx(S, 1)
        mask = hp**2.0 < lmd/beta
        hp[mask] = 0.0

        # with hp, solve for S in Eq. (8)
        S = np.real(ifft((F_I + beta*fft(circulant_dx(hp, -1))) / (1.0 + beta*F_denom)))

        # iteration step
        if return_history:
            S_history.append(np.array(S))
        beta *= beta_rate
        if beta > beta_max: break

    if return_history:
        return S_history

    return S

    def dct(self, x):
        '''Compute discrete cosine transform of 1-d array x'''

        #probably don't need this here anymore since it is in fftpack now    

        N = len(x)
        #calculate DCT weights
        w = (np.exp(-1j*np.arange(N)*np.pi/(2*N))/np.sqrt(2*N))
        w[0] = w[0]/np.sqrt(2)

        #test for odd or even function 
        if (N%2==1) or (any(np.isreal(x)) == False):
            y = np.zeros(2*N)
            y[0:N] = x
            y[N:2*N] = x[::-1]
            yy = fftpack.fft(y)
            yy = yy[0:N]
        else:
            y = np.r_[(x[0:N:2], x[N:0:-2])]
            yy = fftpack.fft(y)
            w = 2*w

        #apply weights
        X = w * yy

        if all(np.isreal(x)):
            X = np.real(X)

        return X

    def compute(self):
        sig_array = self.get_input("Signals")

        # If there is no input on the samples port,
        # use the number of samples in an array row for
        # the number of fft points.
        if self.has_input("Samples"):
            pts = self.get_input("Samples")
            
        else:
            try:
                pts = sig_array.get_shape()[1]
            except:
                pts = sig_array.get_shape()[0]

        sh = sig_array.get_shape()
        if len(sh) < 2:
            shp = (1, sh[0])
            sig_array.reshape(shp)

        (num_sigs, num_samps) = sig_array.get_shape()
        phasors = fftpack.fft(sig_array.get_row_range(0,0), pts)
        out_ar = phasors

        for i in xrange(1,num_sigs):
            phasors = fftpack.fft(sig_array.get_row_range(i,i), pts)
            out_ar = numpy.vstack([out_ar, phasors])
        
        out = NDArray()
        out.set_array(out_ar)
        self.set_output("FFT Output", out)

def fourier_transform_and_reconstruct(image, detrend=False, window=False,
                                      ffunc=None):
    """
    Take fourier transform, alter it, and reconstruct image.  For some
    reason this is shifting the origin by 1 pixel after reconstruction, which
    should not happen.

    :param image: data
    :type image: :py:class:`numpy.ndarray`

    :param ffunc: function that alters FFT matrix
    :type ffunc: func

    :return: modified image data
    :rtype: :py:class:`numpy.ndarray`
    """
    if window:
        w = signal.hamming(image.shape)
    else:
        w = np.ones_like(image)

    if detrend:
        f = fftpack.fft(w * signal.detrend(image))
    else:
        f = fftpack.fft(w * image)

    # alter the fft
    if not ffunc is None:
        f = ffunc(f)

    result = np.fliplr(fftpack.fft(f))

    return result > result.mean()

def xcorrf(trace1, trace2, shift=None):
    """
    Cross-correlation of numpy arrays data1 and data2 in frequency domain.
    """
    data1 = trace1.data
    data2 = trace2.data

    complex_result = data1.dtype == complex or data2.dtype == complex
    N1 = len(data1)
    N2 = len(data2)

    data1 = data1.astype("float64")
    data2 = data2.astype("float64")

    # Always use 2**n-sized FFT, perform xcorr
    size = max(2 * shift + 1, (N1 + N2) // 2 + shift)
    nfft = nextpow2(size)
    IN1 = fft(data1, nfft)
    IN1 *= conjugate(fft(data2, nfft))
    ret = ifft(IN1)
    del IN1

    if not complex_result:
        ret = ret.real
    # shift data for time lag 0 to index 'shift'

    ret = roll(ret, -(N1 - N2) // 2 + shift)[: 2 * shift + 1]

    return copy(ret)

def test_phase_randomize():
    from brainiak.utils.utils import phase_randomize
    import numpy as np
    from scipy.fftpack import fft
    import math
    from scipy.stats import pearsonr

    # Generate auto-correlated signals
    nv = 2
    T = 100
    ns = 3
    D = np.zeros((nv, T, ns))
    for v in range(nv):
        for s in range(ns):
            D[v, :, s] = np.sin(np.linspace(0, math.pi * 5 * (v + 1), T)) + \
                         np.sin(np.linspace(0, math.pi * 6 * (s + 1), T))

    freq = fft(D, axis=1)
    D_pr = phase_randomize(D)
    freq_pr = fft(D_pr, axis=1)
    p_corr = pearsonr(np.angle(freq).flatten(), np.angle(freq_pr).flatten())[0]

    assert np.isclose(abs(freq), abs(freq_pr)).all(), \
        "Amplitude spectrum not preserved under phase randomization"

    assert abs(p_corr) < 0.03, \
        "Phases still correlated after randomization"

def pulseSpectrum(t, SVEAAmp, lambdaZero = 0.0, units = 'nm'):
	'''
	Compute the spectrum of o SVEA pulse center at lambdaZero
		* t: time vector
		* SVEAAmp: SVEA enveloppe of the pulse
		* lambdaZero: center of the pulse [m]
		* units: Units of the output ['nm','um','m']
	'''

	C = 2.99792458e-4
	nt = len(t)
	dt = t[1] - t[0]
	T =  t.max()-t.min()
	w = wspace(T,nt)
	vs = fftpack.fftshift(w/(2*pi))

	# Assign uniScale 
	unitScale = {
	  'nm': lambda: 1.0e9,
	  'um': lambda: 1.0e6,
	  'm': lambda: 1.0
	}[units]()

	if lambdaZero != 0.0:
		wavelength = ( 1.0/( (vs/C)+1.0/(lambdaZero) ) )*unitScale
		return [wavelength, fftpack.fftshift(pow(abs(dt*fftpack.fft(SVEAAmp)/sqrt(2.0*pi)),2))]
	else:
		return [vs, fftpack.fftshift(pow(abs(dt*fftpack.fft(SVEAAmp)/sqrt(2.0*pi)),2))]

def CQT_fast(x,fs,bins,fmin,fmax,M):

	threshold = 0.0054 #for Hamming window
	K = int(bins*np.ceil(np.log2(fmax/fmin)))
	Q = 1/(2**(1/bins)-1)
	nfft = np.int32(nearestPow2(np.ceil(Q*fs/fmin)))
	tempKernel = np.zeros(nfft, dtype = np.complex)
	specKernel = np.zeros(nfft, dtype = np.complex)
	sparKernel = []

	#create sparse Kernel 
	for k in range(K-1,-1,-1):
		fk = (2**(k/bins))*fmin
		N = np.int32(np.round((Q*fs)/fk))
		tempKernel[:N] = hamming(N)/N * np.exp(-2*np.pi*1j*Q*np.arange(N)/N)
		specKernel = fft(tempKernel)
		specKernel[np.where(np.abs(specKernel) <= threshold)] = 0
		if k == K-1:
			sparKernel = specKernel
		else:
			sparKernel = np.vstack((specKernel, sparKernel))
	
	sparKernel = np.transpose(np.conjugate(sparKernel))/nfft
	ft = fft(x,nfft)
	cqt = np.dot(ft, sparKernel)
	ft = fft(x,nfft*(2**M))
	#calculate harmonic power spectrum
	#harm_pow = HPS(ft,M)
	#cqt = np.dot(harm_pow, sparKernel)
	return cqt

    def testAngularMethod(self):
        sample_rate = 16000
        sample_delay = 5
        angle = math.pi / 6
        if abs(math.cos(angle)) > 1e-10:
            dist = sample_delay * pa_tools.SPEED_OF_SOUND / (sample_rate * math.cos(angle))
        else:
            dist = 1
            sample_delay = 0
        print "distance: " + str(dist)
        mics = np.array([[0., 0.], [dist, 0.]], dtype=np.float32)
        data_len = 100
        data1 = np.random.rand(1, data_len)
        if sample_delay > 0:
            data2 = np.concatenate((np.random.rand(1, sample_delay),
                                    [data1[0, :-sample_delay]]), axis=1)
        else:
            data2 = data1
        # Get dfts
        fft1 = fftp.fft(data1[0])
        fft2 = fftp.fft(data2[0])
        ffts = np.array([fft1, fft2])
        loc = DirectionLocalizer(mics, sample_rate=sample_rate)
        direction = loc.get_direction_np(ffts)
        print "direction: " + str(direction)

        # Plot
        plt.figure()
        plt.plot(mics[:, 0], mics[:, 1], 'bo')
        plt.quiver(0, 0, direction[0], direction[1], scale=20)
        plt.show()

    def fft(self, normalize = False):
        """
        params:
            `normalize`: if True, the data will be normalize before pass to fft. Default is False.
        """
        
        coefs = []
        data = self.audio_data.data

        if normalize:
            if self.audio_data.dtype in [_np.uint8, _np.uint16, _np.uint32]:
                data = 2.*data/2**(8*self.audio_data.BIT_WIDTH) - 1
            elif self.audio_data.dtype in [_np.int8, _np.int16, _np.int32]:
                data = 2.*data/2**(4*self.audio_data.BIT_WIDTH) - 1
            else:
                print "[Warning] Unrecognized dtype detected."

        if self.audio_data.CHANNELS == 1: # mono audio data
            # Take only the coef for the positive freq since the data is real-valued.
            N = len(self.audio_data.data)
            c1 = _fftpack.fft(data)[:N/2]
            coefs.append(c1)
        else: # stereo audio data
            N = len(self.audio_data.data[0])
            c1 = _fftpack.fft(data[0])[:N/2]
            coefs.append(c1)
            c2 = _fftpack.fft(data[1])[:N/2]
            coefs.append(c2)
        self.__cached_fft = coefs

        return coefs

def get_phase_diff(target_freq, fs, a, b):
    # sample rate (Hz)
    sample_rate = fs

    # Highest frequency captured - limited by sample rate (Hz)
    high_freq_bound = sample_rate / 2

    # Number samples (bins)
    sample_number = len(a)

    # Signal freq held in each index of arrays (Hz/bin)
    scale = sample_rate / sample_number

    # Index of arrays for complex number we want (bin)
    index = target_freq / scale

    # Update target frequency to one that matches scale of fft
    target_freq = scale * index

    # Getting our phases (radians)
    a_phase = phase(fft(a)[index])
    b_phase = phase(fft(b)[index])
    b_phase = relative_wraparound(a_phase, b_phase)
    print("a_phase = %f pi radians" % (a_phase / math.pi))
    print("b_phase = %f pi radians" % (b_phase / math.pi))

    # Compute phase difference
    phase_diff = b_phase - a_phase
    print("b leads a by %fpi radians" % (phase_diff / math.pi))

    return phase_diff