Python fftpack.ifft函数代码示例

def xcorrnorm(tr1, tr2, pad=True):
    """
    Compute normalized cross correlation of two traces
    maxcor, maxlag, maxdt, cc, lags, tlags = xcorr1x1(tr1, tr2)

    INPUTS
    tr1 - obspy trace1
    tr2 - obspy trace2
    NOT IMPLEMENTED YET
    freqmin, freqmax - optional, restrict frequency range to [freqmin, freqmax]
    lags = lag to compute if None, will compute all lags and find max,

    OUTPUTS
    maxcor - value of maximum correlation
    maxlag - lag of maximum correlation (in samples) - this is the number of samples to shift tr2 so it lines up with tr1
    maxdt - time lag of max lag, in seconds
    cc - cross correlation value at each shift
    lags - lag at each cc value in samples
    tlags - lag at each cc value in seconds

    TODO
    add option to only compute certain lags
    add option to output entire cc function
    """
    from scipy.fftpack import fft, ifft
    if tr1.stats.sampling_rate != tr2.stats.sampling_rate:
        raise RuntimeError('tr1 and tr2 have different sampling rates')
        return
    # make sure data is float
    dat1 = tr1.data*1.
    dat2 = tr2.data*1.
    if len(tr1) != len(tr2):
        if pad is True:
            print('tr1 and tr2 have different lengths, padding with zeros')
            if len(dat1) > len(dat2):
                dat2 = np.lib.pad(dat2, (0, len(dat1)-len(dat2)), 'constant', constant_values=(0., 0.))
            else:
                dat1 = np.lib.pad(dat1, (0, len(dat2)-len(dat1)), 'constant', constant_values=(0., 0.))
        else:
            raise RuntimeError('tr1 and tr2 are different lengths, set pad=True if you want to proceed')
            return
    # pad data to double number of samples to avoid wrap around and pad more to next closest power of 2 for fft
    n2 = nextpow2(len(dat1))

    FFT1 = fft(dat1, n2)
    norm1 = np.sqrt(np.real(ifft(FFT1*np.conj(FFT1), n2)))
    FFT2 = fft(dat2, n2)
    norm2 = np.sqrt(np.real(ifft(FFT2*np.conj(FFT2), n2)))
    cctemp = np.real(ifft(FFT1*np.conj(FFT2), n2))
    cc = cctemp/(norm1[0]*norm2[0])
    M = len(FFT1)
    lags = np.roll(np.linspace(-M/2 + 1, M/2, M, endpoint=True), M/2 + 1).astype(int)
    indx = np.argmax(cc)

    maxcor = cc[indx]
    maxlag = lags[indx]
    maxdt = 1./tr1.stats.sampling_rate*maxlag
    tlags = 1./tr1.stats.sampling_rate*lags

    return maxcor, maxlag, maxdt, cc, lags, tlags

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 ifcglt(A): # Lobatto Nodal to Modal Coefficients
    """
    Fast Chebyshev-Gauss-Lobatto transformation from 
    point space values (nodal) to Chebyshev expansion 
    coefficients (modal). If I=numpy.identity(n), then
    Ti=chebyshev.ifcglt(I) will be the inverse of the 
    Chebyshev Vandermonde matrix on the Lobatto nodes
    """
    size = A.shape
    m = size[0]
    k = m-1-np.arange(m-1)
  
    if len(size) == 2: # Multiple vectors
        V = np.vstack((A[0:m-1,:],A[k,:]))
        F = ifft(V, n=None, axis=0)
        B = np.vstack((F[0,:],2*F[1:m-1,:],F[m-1,:]))
    else:  # Single vector
        V = np.hstack((A[0:m-1],A[k]))
        F = ifft(V, n=None)
        B = np.hstack((F[0],2*F[1:m-1],F[m-1]))
        
    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 XIntegralsFFT(GF_A,Bubble_A,Lambda,BubZero):
	''' calculate X integral to susceptibilities using FFT '''
	N = int((len(En_A)-1)/2)
	Kappa_A  = TwoParticleBubble(GF_A,GF_A**2,'eh')
	Bubble_A = TwoParticleBubble(GF_A,GF_A,'eh')
	#print(Kappa_A[N],Bubble_A[N])
	V_A   = 1.0/(1.0+Lambda*Bubble_A)
	KV_A  = Lambda*Kappa_A*V_A**2
	KmV_A = Lambda*sp.flipud(sp.conj(Kappa_A))*V_A**2
	## zero-padding the arrays
	exFD_A  = sp.concatenate([FD_A[N:],sp.zeros(2*N+2),FD_A[:N+1]])
	ImGF_A  = sp.concatenate([sp.imag(GF_A[N:]),sp.zeros(2*N+2),sp.imag(GF_A[:N+1])])
	ImGF2_A = sp.concatenate([sp.imag(GF_A[N:]**2),sp.zeros(2*N+2),sp.imag(GF_A[:N+1]**2)])
	ImV_A   = sp.concatenate([sp.imag(V_A[N:]),sp.zeros(2*N+2),sp.imag(V_A[:N+1])])
	ImKV_A  = sp.concatenate([sp.imag(KV_A[N:]),sp.zeros(2*N+2),sp.imag(KV_A[:N+1])])
	ImKmV_A = sp.concatenate([sp.imag(KmV_A[N:]),sp.zeros(2*N+2),sp.imag(KmV_A[:N+1])])
	## performing the convolution
	ftImX11_A = -sp.conj(fft(exFD_A*ImV_A))*fft(ImGF2_A)*dE
	ftImX12_A =  fft(exFD_A*ImGF2_A)*sp.conj(fft(ImV_A))*dE
	ftImX21_A = -sp.conj(fft(exFD_A*ImKV_A))*fft(ImGF_A)*dE
	ftImX22_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKV_A))*dE
	ftImX31_A = -sp.conj(fft(exFD_A*ImKmV_A))*fft(ImGF_A)*dE
	ftImX32_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKmV_A))*dE
	## inverse transform
	ImX1_A =  sp.real(ifft(ftImX11_A+ftImX12_A))/sp.pi
	ImX2_A =  sp.real(ifft(ftImX21_A+ftImX22_A))/sp.pi
	ImX3_A = -sp.real(ifft(ftImX31_A+ftImX32_A))/sp.pi
	ImX1_A =  sp.concatenate([ImX1_A[3*N+4:],ImX1_A[:N+1]])
	ImX2_A =  sp.concatenate([ImX2_A[3*N+4:],ImX2_A[:N+1]])
	ImX3_A =  sp.concatenate([ImX3_A[3*N+4:],ImX3_A[:N+1]])
	## getting real part from imaginary
	X1_A = KramersKronigFFT(ImX1_A) + 1.0j*ImX1_A + BubZero # constant part !!!
	X2_A = KramersKronigFFT(ImX2_A) + 1.0j*ImX2_A
	X3_A = KramersKronigFFT(ImX3_A) + 1.0j*ImX3_A
	return [X1_A,X2_A,X3_A]

    def bench_random(self):
        from numpy.fft import ifft as numpy_ifft
        print()
        print('       Inverse Fast Fourier Transform')
        print('===============================================')
        print('      |     real input    |    complex input   ')
        print('-----------------------------------------------')
        print(' size |  scipy  |  numpy  |  scipy  |  numpy  ')
        print('-----------------------------------------------')
        for size,repeat in [(100,7000),(1000,2000),
                            (256,10000),
                            (512,10000),
                            (1024,1000),
                            (2048,1000),
                            (2048*2,500),
                            (2048*4,500),
                            ]:
            print('%5s' % size, end=' ')
            sys.stdout.flush()

            for x in [random([size]).astype(double),
                      random([size]).astype(cdouble)+random([size]).astype(cdouble)*1j
                      ]:
                if size > 500: y = ifft(x)
                else: y = direct_idft(x)
                assert_array_almost_equal(ifft(x),y)
                print('|%8.2f' % measure('ifft(x)',repeat), end=' ')
                sys.stdout.flush()

                assert_array_almost_equal(numpy_ifft(x),y)
                print('|%8.2f' % measure('numpy_ifft(x)',repeat), end=' ')
                sys.stdout.flush()

            print(' (secs for %s calls)' % (repeat))
        sys.stdout.flush()

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 idct(self, X):
        '''Compute inverse discrete cosine transform of a 1-d array X'''
        
        N = len(X)
        w = np.sqrt(2*N)*np.exp(1j*np.arange(N)*np.pi/(2.*N))

        if (N%2==1) or (any(np.isreal(X))== False):
            w[0] = w[0] * np.sqrt(2)
            yy = np.zeros(2*N)
            yy[0:N] = w * X
            yy[N+1:2*N] = -1j * w[1:N] * X[1:N][::-1]
            y = fftpack.ifft(yy)
            x = y[0:N]
        else:
            w[0] = w[0]/np.sqrt(2)
            yy = X *w
            y = fftpack.ifft(yy)
            x = np.zeros(N)
            x[0:N:2] = y[0:N/2]
            x[1:N:2] = y[N:(N/2)-1:-1]

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

        return x

def convolve_power(power_spectrum, window_power, axis=-1):
    """Convolve a power spectrum with a window function."""
    
    data_corr = fft.ifft(power_spectrum, axis=axis)
    window_corr = fft.ifft(window_power, axis=axis)
    true_corr = data_corr * window_corr
    true_power = fft.fft(true_corr, axis=axis, overwrite_x=True)
    return true_power

def sprModel(x, fs, w, N, t):
	"""
	Analysis/synthesis of a sound using the sinusoidal plus residual model, one frame at a time
	x: input sound, fs: sampling rate, w: analysis window, 
	N: FFT size (minimum 512), t: threshold in negative dB, 
	returns y: output sound, ys: sinusoidal component, xr: residual component
	"""

	hN = N/2                                                      # size of positive spectrum
	hM1 = int(math.floor((w.size+1)/2))                           # half analysis window size by rounding
	hM2 = int(math.floor(w.size/2))                               # half analysis window size by floor
	Ns = 512                                                      # FFT size for synthesis (even)
	H = Ns/4                                                      # Hop size used for analysis and synthesis
	hNs = Ns/2      
	pin = max(hNs, hM1)                                           # initialize sound pointer in middle of analysis window          
	pend = x.size - max(hNs, hM1)                                 # last sample to start a frame
	fftbuffer = np.zeros(N)                                       # initialize buffer for FFT
	ysw = np.zeros(Ns)                                            # initialize output sound frame
	xrw = np.zeros(Ns)                                            # initialize output sound frame
	ys = np.zeros(x.size)                                         # initialize output array
	xr = np.zeros(x.size)                                         # initialize output array
	w = w / sum(w)                                                # normalize analysis window
	sw = np.zeros(Ns)     
	ow = triang(2*H)                                              # overlapping window
	sw[hNs-H:hNs+H] = ow      
	bh = blackmanharris(Ns)                                       # synthesis window
	bh = bh / sum(bh)                                             # normalize synthesis window
	wr = bh                                                       # window for residual
	sw[hNs-H:hNs+H] = sw[hNs-H:hNs+H] / bh[hNs-H:hNs+H]
	while pin<pend:  
  #-----analysis-----             
		x1 = x[pin-hM1:pin+hM2]                                     # select frame
		mX, pX = DFT.dftAnal(x1, w, N)                              # compute dft
		ploc = UF.peakDetection(mX, t)                              # find peaks 
		iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)         # refine peak values		iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)          # refine peak values
		ipfreq = fs*iploc/float(N)                                  # convert peak locations to Hertz
		ri = pin-hNs-1                                              # input sound pointer for residual analysis
		xw2 = x[ri:ri+Ns]*wr                                        # window the input sound                                       
		fftbuffer = np.zeros(Ns)                                    # reset buffer
		fftbuffer[:hNs] = xw2[hNs:]                                 # zero-phase window in fftbuffer
		fftbuffer[hNs:] = xw2[:hNs]                           
		X2 = fft(fftbuffer)                                         # compute FFT for residual analysis
  #-----synthesis-----
		Ys = UF.genSpecSines(ipfreq, ipmag, ipphase, Ns, fs)        # generate spec of sinusoidal component          
		Xr = X2-Ys;                                                 # get the residual complex spectrum
		fftbuffer = np.zeros(Ns)
		fftbuffer = np.real(ifft(Ys))                               # inverse FFT of sinusoidal spectrum
		ysw[:hNs-1] = fftbuffer[hNs+1:]                             # undo zero-phase window
		ysw[hNs-1:] = fftbuffer[:hNs+1] 
		fftbuffer = np.zeros(Ns)
		fftbuffer = np.real(ifft(Xr))                               # inverse FFT of residual spectrum
		xrw[:hNs-1] = fftbuffer[hNs+1:]                             # undo zero-phase window
		xrw[hNs-1:] = fftbuffer[:hNs+1]
		ys[ri:ri+Ns] += sw*ysw                                      # overlap-add for sines
		xr[ri:ri+Ns] += sw*xrw                                      # overlap-add for residual
		pin += H                                                    # advance sound pointer
	y = ys+xr                                                     # sum of sinusoidal and residual components
	return y, ys, xr

def ImpulseResponse(H, F):
    """
    Input:
        H = existing transferfunction spectrum
        F = [Hz] Frequency bins array
    output:
        IR = [sec] time signal of resultin impulse respose
        T = time of the total impusle response inluding silence
        fs = measured sample frequency

    Impulse response makes use of ifft method to calculate a time signal back
    from the transferfunction.
    fft of the form:
               n-1                      m*k
    a_m =1/n * sum A_k * exp ( -2pi * i ----)    for m = 0, ..., n-1
               k=0                       n
    Therefore fft * 1/N to correct amplitude

    TODO:
    - correction for N
    - Silence removal
    """
    print("""This function works only correct when:
            - no smoothing or averaging is applied
            - full spectrum data is returned!""")
    from scipy.fftpack import ifft, fftshift, fftfreq
    # 2DO Check for full Spectrum
    if isinstance(H, tuple):  # could also use assert
        # Check arrays for type of signal
        # Mag + Phase --> Mag is absolut no neg values
        # Check symetrie no symmerie is false info...
        # Without shift X[0] = F - 0Hz and X[fs/2] =
        # shift check on symmetry?
        H0even = Symmetry(H[0], 'even')
        H1odd = Symmetry(H[1], 'odd')
        if (H0even and H1odd):
            if H[0] == abs(H[0]):
                if max(abs(H[1])) <= np.pi:  # check for phase RAD information
                    pass
                elif max(abs(H[1])) <= 180:  # check for phase DEG information
                    H[1] = H[1] * np.pi / 180
                    pass
#                    http://stackoverflow.com/questions/2598734/numpy-creating-a-complex-array-from-2-real-ones
#            elif ...:  # check for complex paterns don't know how
#                pass
            else:
                pass
        else:
            raise MeasError.DataError(H, 'No valid frequency array')
    elif np.any(np.iscomplex(H)):
        IR = ifft(H)
    else:
        raise TypeError('Not right input type')
    IR = ifft(H)
    dF = F[1]-F[0]
    fs = np.round(len(F) * dF)
    T = len(F) * 1 / fs
    return(IR, fs, T)

    def test_definition(self):
        x = np.array([1, 2, 3, 4 + 1j, 1, 2, 3, 4 + 2j], self.cdt)
        y = ifft(x)
        y1 = direct_idft(x)
        assert_equal(y.dtype, self.cdt)
        assert_array_almost_equal(y, y1)

        x = np.array([1, 2, 3, 4 + 0j, 5], self.cdt)
        assert_array_almost_equal(ifft(x), direct_idft(x))

 def test_random_real(self):
     for size in [1, 51, 111, 100, 200, 64, 128, 256, 1024]:
         x = random([size]).astype(self.rdt)
         y1 = ifft(fft(x))
         y2 = fft(ifft(x))
         assert_equal(y1.dtype, self.cdt)
         assert_equal(y2.dtype, self.cdt)
         assert_array_almost_equal(y1, x)
         assert_array_almost_equal(y2, x)

def rceps(x): 
	y = sp.real(ifft(sp.log(sp.absolute(fft(x)))))
	n = len(x) 
	if (n%2) == 1:
		ym = np.hstack((y[0], 2*y[1:n/2], np.zeros(n/2-1)))
	else:
		ym = np.hstack((y[0], 2*y[1:n/2], y[n/2+1], np.zeros(n/2-1)))
	ym = sp.real(ifft(sp.exp(fft(ym)))) 
	return (y, ym)

def ACF_scargle(t, y, dy, n_omega=2 ** 10, omega_max=100):
    """Compute the Auto-correlation function via Scargle's method

    Parameters
    ----------
    t : array_like
        times of observation.  Assumed to be in increasing order.
    y : array_like
        values of each observation.  Should be same shape as t
    dy : float or array_like
        errors in each observation.
    n_omega : int (optional)
        number of angular frequencies at which to evaluate the periodogram
        default is 2^10
    omega_max : float (optional)
        maximum value of omega at which to evaluate the periodogram
        default is 100

    Returns
    -------
    ACF, t : ndarrays
        The auto-correlation function and associated times
    """
    t = np.asarray(t)
    y = np.asarray(y)

    if y.shape != t.shape:
        raise ValueError("shapes of t and y must match")

    dy = np.asarray(dy) * np.ones(y.shape)

    d_omega = omega_max * 1. / (n_omega + 1)
    omega = d_omega * np.arange(1, n_omega + 1)

    # recall that P(omega = 0) = (chi^2(0) - chi^2(0)) / chi^2(0)
    #                          = 0

    # compute P and shifted full-frequency array
    P = lomb_scargle(t, y, dy, omega,
                     generalized=True)
    P = np.concatenate([[0], P, P[-2::-1]])

    # compute PW, the power of the window function
    PW = lomb_scargle(t, np.ones(len(t)), dy, omega,
                      generalized=False, subtract_mean=False)
    PW = np.concatenate([[0], PW, PW[-2::-1]])

    # compute the  inverse fourier transform of P and PW
    rho = fftpack.ifft(P).real
    rhoW = fftpack.ifft(PW).real

    ACF = fftpack.fftshift(rho / rhoW) / np.sqrt(2)
    N = len(ACF)
    dt = 2 * np.pi / N / (omega[1] - omega[0])
    t = dt * (np.arange(N) - N // 2)

    return ACF, t