Python fftpack.fftshift函数代码示例

def rescale_target_superpixel_resolution(E_target):
    '''Rescale the target field to the superpixel resolution (currently only 4x4 superpixels implemented)'''
    
    superpixelSize = 4
    
    ny,nx = scipy.shape(E_target)
    
    maskCenterX = scipy.ceil((nx+1)/2)
    maskCenterY = scipy.ceil((ny+1)/2)
    
    nSuperpixelX = int(nx/superpixelSize)
    nSuperpixelY = int(ny/superpixelSize)
    
    FourierMaskSuperpixelResolution = fourier_mask(ny,nx,superpixelSize)
    
    
    E_target_ft = fft.fftshift(fft.fft2(fft.ifftshift(E_target)))
    
    #Apply mask
    E_target_ft = FourierMaskSuperpixelResolution*E_target_ft
    
    #Remove zeros outside of mask
    E_superpixelResolution_ft = E_target_ft[(maskCenterY - scipy.ceil((nSuperpixelY-1)/2)-1):(maskCenterY + scipy.floor((nSuperpixelY-1)/2)),(maskCenterX - scipy.ceil((nSuperpixelX-1)/2)-1):(maskCenterX + scipy.floor((nSuperpixelX-1)/2))]
    
    # Add phase gradient to compensate for anomalous 1.5 pixel shift in real
    # plane
    
    phaseFactor = [[(scipy.exp(2*1j*pi*((k+1)/nSuperpixelY+(j+1)/nSuperpixelX)*3/8)) for j in range(nSuperpixelX)] for k in range(nSuperpixelY)] # QUESTION
    E_superpixelResolution_ft = E_superpixelResolution_ft*phaseFactor
    
    # Fourier transform back to DMD plane
    E_superpixelResolution = fft.fftshift(fft.ifft2(fft.ifftshift(E_superpixelResolution_ft)))

    return E_superpixelResolution

 def step4(self):
   '''
   Perform a 4th order timestep
   '''
   
   def order2(c):
     Vc = np.exp( -1j * c * self.dt / 2. * 
                   ( self.V - self.gravity() + 
                     self.g * abs( self.psi ) ** 2
                   )
                 )
     Tc = self.expksquare ** c
     return Vc, Tc
   
   p = 1/(4.-4.**(1/3.))
   q = 1 - 4 * p
   
   Vp,Tp = order2( p )
   Vq,Tq = order2( q )
   
   return Vp * ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp ** 2 * 
               ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp * Vq *
               ff.fftshift( ff.ifft2( Tq * ff.fft2( ff.fftshift( Vq * Vp * 
               ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp ** 2 *  
               ff.fftshift( ff.ifft2( Tp * ff.fft2( ff.fftshift( Vp * self.psi 
               ) ) ) )
               ) ) ) )
               ) ) ) )
               ) ) ) )
               ) ) ) )

def compute_fft(x, n=None, axis=True, power=True, fftshift=True, fs=1.0):
    """
    Compute and return the FFT of the input data
    :param x:
    :param n:
    :param power:
    :param fftshift:
    :return:
    """
    x = x.strip('[]')
    x = [float(i) for i in x.split(',') if i is not '']

    fft_out = fftpack.fft(x,n)

    if axis is True:
        axis_values = np.linspace(0, 2*fs, len(fft_out))
    else:
        axis_values = []

    if power is True:
        fft_out = np.abs(fft_out)

    if fftshift is True:
        fft_out = fftpack.fftshift(fft_out)
        axis_values = fftpack.fftshift(axis_values)

    outdict = {'data_in': x,
               'data_out': fft_out,
               'axis':axis_values,
                'power':power,
                'fft_shift':fft_shift}

    def plot_transfer_function(self, ds):
        ax1 = pl.subplot(211)
        ax2 = pl.subplot(212)

        Hw = self.transfer_function(ds.omega)
        ax1.plot(fftp.fftshift(ds.f), np.abs(fftp.fftshift(Hw)) ** 2, "k--")
        ax2.plot(fftp.fftshift(ds.f), np.angle(fftp.fftshift(Hw)), "k--")

def test5():
    global L0, N

    L = deepcopy(L0)
    rho = zeros(N, 'double')
    rho[0]  = 1.
    rho[N/2] = 1.

    print rho

    print fft(rho)
    rho = fftshift(fft(rho))
    print "fft(rho) =", rho

    L = fft(L).T
    L = fft(L).T

    L = fftshift(L)

    #print L

    x = linalg.solve(L, rho)
    print "x =", x
    #x[abs(x)<0.001] = 0

    x = ifftshift(ifft(x)).real * N
    print "ifft(x) =", x
    F = []
    for i in xrange(len(x)-1):
        F.append(x[i+1] - x[i])

    print "F =", F
    print "--------------------------------"

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 rolloff(grid_shape, nhood_shape=[6,6], ncrop=None, osfactor=2, threshold=0.01, axes=[-1, -2], combi=False, slice_prof_coef=1):
    if combi:
        nz = spatial_dims[-1]
        spatial_dims = grid_shape[:-1]
        ndims = len(spatial_dims)
        ro_data = ones(nz, dtype='complex64')
        ro_K = zeros((nz, ndims), dtype='float32')
        ro_K[:,2] = linspace(-nz/2, nz/2, nz, endpoint=False)
    else:
        spatial_dims = grid_shape
        ndims = len(spatial_dims)
        ro_data = array([1.0], dtype='complex64')
        ro_K = array([[0]*ndims], dtype='float32')
    G = sparse_matrix_operator(ro_K, grid_shape=spatial_dims, nhood_shape=nhood_shape, osfactor=osfactor, combi=combi, slice_prof_coef=slice_prof_coef)
    ro = G.T * ro_data
    n = sqrt(G.shape[1])
    ro = reshape(ro, (n,n))
    #if osfactor > 1 and ncrop == None:
    #    ncrop = int((spatial_dims[0]/osfactor) * (osfactor - 1) / 2)
    ro = fftshift(abs(ifftn(fftshift(ro, axes=axes), axes=axes)), axes=axes)  # transform to image
    if ncrop > 0:
        ro = ro[ncrop:-ncrop, ncrop:-ncrop]
    #print 'rolloff shape:', ro.shape
    ro = ro / ro.max()  # normalize
    ro[ro < threshold] = 1.0
    #ro_max = ro.max()
    #ro[ro < threshold*ro_max] = 1.0
    ro = 1.0 / ro
    #ro = ro**2
    #print 'TOOK OUT SQUARED RO'
    #ro = ro / ro.max()
    return ro

    def __init__(self, x, y, s, detrend=True, window=False, **kwargs):
        # r-space
        self.x  = np.asanyarray(x)
        self.y  = np.asanyarray(y)
        self.s  = np.asanyarray(s)

        assert len(self.x.shape) == 2
        assert self.x.shape == self.y.shape == self.s.shape
        assert self.x.size == self.y.size == self.s.size

        # r-space spacing
        self.dx = self._delta(self.x, np.index_exp[0,0], np.index_exp[1,0])
        self.dy = self._delta(self.y, np.index_exp[0,0], np.index_exp[0,1])

        # r-space samples
        self.n0 = self.x.shape[0]
        self.n1 = self.x.shape[1]

        # r-space lengths
        self.lx = self.n0 * self.dx
        self.ly = self.n1 * self.dy

        # k-space
        u = fftpack.fftshift(fftpack.fftfreq(self.n0))
        v = fftpack.fftshift(fftpack.fftfreq(self.n1))
        self.u, self.v = np.meshgrid(u, v, indexing='ij')

        # k-space spacing
        self.du = self._delta(self.u, np.index_exp[0,0], np.index_exp[1,0])
        self.dv = self._delta(self.v, np.index_exp[0,0], np.index_exp[0,1])

        # k-space lengths
        self.lu = self.n0 * self.du
        self.lv = self.n1 * self.dv

        # nyquist
        try:
            self.nyquist_u = 0.5/self.dx
        except ZeroDivisionError:
            self.nyquist_u = 0.0

        try:
            self.nyquist_v = 0.5/self.dy
        except ZeroDivisionError:
            self.nyquist_v = 0.0

        self.k = np.sqrt(self.u**2 + self.v**2)

        # detrend the signal
        if detrend:
            self.s = signal.detrend(self.s)

        # apply window to signal
        if window:
            self._window()
            self.s = self.s * self.window

        # compute the FFT
        self.fft = fftpack.fftshift(fftpack.fft2(self.s))
        self.power = self.fft.real**2 + self.fft.imag**2

def genwavenumber(nlon):
    if (nlon%2 == 0):
        wavenumber = fftpack.fftshift(fftpack.fftfreq(nlon)*nlon)[1:]
    else:
        wavenumber = fftpack.fftshift(fftpack.fftfreq(nlon)*nlon)

    return wavenumber

    def fft(self, nfft=None, ssb=False):
        """
        Computes the Fast Fourier transform of the signal using :func:`scipy.fftpack.fft` function.
        The Fourier transform of a time series function is defined as:

        .. math::
           \mathcal{F}(y) ~=~ \int_{-\infty}^{\infty} y(t) e^{-2 \pi j f t}\,dt

        Parameters
        ----------
        nfft : int, optional
            Specify the number of points for the FFT. The default is the length of
            the time series signal.

        ssb : boolean, optional
            If true, returns only the single side band (components corresponding to positive
            frequency).

        Returns
        -------
        : Signal
            The FFT of the signal.
        """
        if nfft is None:
            nfft = self.size
        uf = Signal(fftshift(fft(self.values, n=nfft)), index=fftshift(fftfreq(nfft, self.ts)))
        return uf[uf.index >= 0] if ssb else uf

def sineSubtraction(x, N, H, sfreq, smag, sphase, fs):
	"""
	Subtract sinusoids from a sound
	x: input sound, N: fft-size, H: hop-size
	sfreq, smag, sphase: sinusoidal frequencies, magnitudes and phases
	returns xr: residual sound
	"""

	hN = N/2                                           # half of fft size
	x = np.append(np.zeros(hN),x)                      # add zeros at beginning to center first window at sample 0
	x = np.append(x,np.zeros(hN))                      # add zeros at the end to analyze last sample
	bh = blackmanharris(N)                             # blackman harris window
	w = bh/ sum(bh)                                    # normalize window
	sw = np.zeros(N)                                   # initialize synthesis window
	sw[hN-H:hN+H] = triang(2*H) / w[hN-H:hN+H]         # synthesis window
	L = sfreq.shape[0]                                 # number of frames, this works if no sines
	xr = np.zeros(x.size)                              # initialize output array
	pin = 0
	for l in range(L):
		xw = x[pin:pin+N]*w                              # window the input sound
		X = fft(fftshift(xw))                            # compute FFT
		Yh = UF_C.genSpecSines(N*sfreq[l,:]/fs, smag[l,:], sphase[l,:], N)   # generate spec sines
		Xr = X-Yh                                        # subtract sines from original spectrum
		xrw = np.real(fftshift(ifft(Xr)))                # inverse FFT
		xr[pin:pin+N] += xrw*sw                          # overlap-add
		pin += H                                         # advance sound pointer
	xr = np.delete(xr, range(hN))                      # delete half of first window which was added in stftAnal
	xr = np.delete(xr, range(xr.size-hN, xr.size))     # delete half of last window which was added in stftAnal
	return xr

 def sim_pic(self,data,alpha):
     '''Do the forward transform to simulate a picture. Currently with 4Pi cruft.'''
     self.alpha = alpha
     self.e1 = fftshift(exp(1j*self.alpha))
     self.e2 = fftshift(exp(2j*self.alpha))
     
     return self.Afunc(data)

def plot_q_qhat(q, t):
    # Plot Potential Vorticity
    plt.clf()
    plt.subplot(2,1,1)
    plt.pcolormesh(xx/1e3,yy/1e3,q)
    plt.colorbar()
    plt.axes([-Lx/2e3, Lx/2e3, -Ly/2e3, Ly/2e3])
    name = "PV at t = %5.2f" % (t/(3600.0*24.0))
    plt.title(name)

    # compute power spectrum and shift ffts
    qe = np.vstack((q0,-np.flipud(q)))
    qhat = np.absolute(fftn(qe))
    kx = fftshift((parms.ikx/parms.ikx[0,1]).real)
    ky = fftshift((parms.iky/parms.iky[1,0]).real)
    qhat = fftshift(qhat)

    Sx, Sy = int(parms.Nx/2), parms.Ny
    Sk = 1.5

    # Plot power spectrum
    plt.subplot(2,1,2)
    #plt.pcolor(kx[Sy:Sy+20,Sx:Sx+20],ky[Sy:Sy+20,Sx:Sx+20],qhat[Sy:Sy+20,Sx:Sx+20])
    plt.pcolor(kx[Sy:int(Sk*Sy),Sx:int(Sk*Sx)],ky[Sy:int(Sk*Sy),Sx:int(Sk*Sx)],
               qhat[Sy:int(Sk*Sy),Sx:int(Sk*Sx)])
    plt.axis([0, 10, 0, 10])
    plt.colorbar()
    name = "PS at t = %5.2f" % (t/(3600.0*24.0))
    plt.title(name)

    plt.draw()

    def test_depth_kwarg(self):
        nz = 100
        zN2 = 0.5*nz**-1 + np.arange(nz, dtype=np.float64)/nz
        zU = 0.5*nz**-1 + np.arange(nz+1, dtype=np.float64)/nz
        N2 = np.full(nz, 1.)
        f0 = 1.
        #beta = 1e-6
        beta = 0.
        Nx = 1
        Ny = 1
        dx = .1
        dy = .1
        k = fft.fftshift( fft.fftfreq(Nx, dx) )
        l = fft.fftshift( fft.fftfreq(Ny, dy) )
        ubar = np.zeros(nz+1)
        vbar = np.zeros(nz+1)
        etax = np.zeros(2)
        etay = np.zeros(2)

        with self.assertRaises(ValueError):
            z, growth_rate, vertical_modes = modes.instability_analysis_from_N2_profile(
                zN2, N2, f0, beta, k, l, zU, ubar, vbar, etax, etay, depth=zN2[-5]
            )

        # no error expected
        z, growth_rate, vertical_mode = modes.instability_analysis_from_N2_profile(
                zN2, N2, f0, beta, k, l, zU, ubar, vbar, etax, etay, depth=1.1
        )

def standard_dsi_algorithm(S,bvals,bvecs):
    #volume size
    sz=16
    #shifting
    origin=8
    #hanning width
    filter_width=32.
    #number of signal sampling points
    n=515

    #odf radius
    #radius=np.arange(2.1,30,.1)
    radius=np.arange(2.1,6,.2)
    #radius=np.arange(.1,6,.1)   
    
    bv=bvals
    bmin=np.sort(bv)[1]
    bv=np.sqrt(bv/bmin)
    qtable=np.vstack((bv,bv,bv)).T*bvecs
    qtable=np.floor(qtable+.5)
   
    #calculate radius for the hanning filter
    r = np.sqrt(qtable[:,0]**2+qtable[:,1]**2+qtable[:,2]**2)
        
    #setting hanning filter width and hanning
    hanning=.5*np.cos(2*np.pi*r/filter_width)
    
    #center and index in q space volume
    q=qtable+origin
    q=q.astype('i8')
    
    #apply the hanning filter
    values=S*hanning
    
    #create the signal volume    
    Sq=np.zeros((sz,sz,sz))
    for i in range(n):        
        Sq[q[i][0],q[i][1],q[i][2]]+=values[i]
    
    #apply fourier transform
    Pr=fftshift(np.abs(np.real(fftn(fftshift(Sq),(sz,sz,sz)))))

    #vertices, edges, faces  = create_unit_sphere(5)    
    #vertices, faces = sphere_vf_from('symmetric362')           
    vertices, faces = sphere_vf_from('symmetric724')           
    odf = np.zeros(len(vertices))
        
    for m in range(len(vertices)):
        
        xi=origin+radius*vertices[m,0]
        yi=origin+radius*vertices[m,1]
        zi=origin+radius*vertices[m,2]
        
        PrI=map_coordinates(Pr,np.vstack((xi,yi,zi)),order=1)
        for i in range(len(radius)):
            odf[m]=odf[m]+PrI[i]*radius[i]**2
   
    peaks,inds=peak_finding(odf.astype('f8'),faces.astype('uint16'))

    return Pr,odf,peaks