Python scipy.fftpack.fftshift() Examples

def from_recip(y):
    """
    Converts Fourier frequencies to spatial coordinates.

    Parameters
    ----------
    y : `list` [`numpy.ndarray` [`float`]], of shape [(nx,), (ny,), ...]
        List (or equivalent) of vectors which define a mesh in the dimension
        equal to the length of `x`

    Returns
    -------
    x : `list` [`numpy.ndarray` [`float`]], of shape [(nx,), (ny,), ...]
        List of vectors defining a mesh such that for a function, `f`, defined on
        the mesh given by `y`, ifft(f) is defined on the mesh given by `x`. 0 will be
        in the middle of `x`.
    """
    x = []
    for Y in y:
        if Y.size > 1:
            x.append(fftfreq(Y.size, Y.item(1) - Y.item(0)) * (2 * pi))
        else:
            x.append(array([0]))
        x[-1] = x[-1].astype(Y.dtype, copy=False)
    return [fftshift(X) for X in x] 

def get_scipy(shape, fftn_shape=None, **kwargs):
    import numpy.fft as numpy_fft
    import scipy.fftpack as scipy_fft

    # use numpy implementation of rfft2/irfft2 because they have not been
    # implemented in scipy.fftpack
    f = {
        "fft2": scipy_fft.fft2,
        "ifft2": scipy_fft.ifft2,
        "rfft2": numpy_fft.rfft2,
        "irfft2": lambda X: numpy_fft.irfft2(X, s=shape),
        "fftshift": scipy_fft.fftshift,
        "ifftshift": scipy_fft.ifftshift,
        "fftfreq": scipy_fft.fftfreq,
    }
    if fftn_shape is not None:
        f["fftn"] = scipy_fft.fftn
    fft = SimpleNamespace(**f)

    return fft 

def get_numpy(shape, fftn_shape=None, **kwargs):
    import numpy.fft as numpy_fft

    f = {
        "fft2": numpy_fft.fft2,
        "ifft2": numpy_fft.ifft2,
        "rfft2": numpy_fft.rfft2,
        "irfft2": lambda X: numpy_fft.irfft2(X, s=shape),
        "fftshift": numpy_fft.fftshift,
        "ifftshift": numpy_fft.ifftshift,
        "fftfreq": numpy_fft.fftfreq,
    }
    if fftn_shape is not None:
        f["fftn"] = numpy_fft.fftn
    fft = SimpleNamespace(**f)

    return fft 

def to_recip(x):
    """
    Converts spatial coordinates to Fourier frequencies.

    Parameters
    ----------
    x : `list` [`numpy.ndarray` [`float`]], of shape [(nx,), (ny,), ...]
        List (or equivalent) of vectors which define a mesh in the dimension
        equal to the length of `x`

    Returns
    -------
    y : `list` [`numpy.ndarray` [`float`]], of shape [(nx,), (ny,), ...]
        List of vectors defining a mesh such that for a function, `f`, defined on
        the mesh given by `x`, `fft(f)` is defined on the mesh given by `y`
    """
    y = []
    for X in x:
        if X.size > 1:
            y.append(fftfreq(X.size, X.item(1) - X.item(0)) * (2 * pi))
        else:
            y.append(array([0]))
        y[-1] = y[-1].astype(X.dtype, copy=False)
    return [fftshift(Y) for Y in y] 

def SO3_ifft(f_hat):
    """
    """
    b = len(f_hat)
    d = setup_d_transform(b)

    df_hat = [d[l] * f_hat[l][:, None, :] for l in range(len(d))]

    # Note: the frequencies where m=-B or n=-B are set to zero,
    # because they are not used in the forward transform either
    # (the forward transform is up to m=-l, l<B
    F = np.zeros((2 * b, 2 * b, 2 * b), dtype=complex)
    for l in range(b):
        F[b - l:b + l + 1, :,  b - l:b + l + 1] += df_hat[l]

    F = fftshift(F, axes=(0, 2))
    f = ifft2(F, axes=(0, 2))
    return f * 2 * (b ** 2) / np.pi 

def deconvolve(f, g):
    """
    FFT based deconvolution

    :param f: array
    :param g: array
    :return: array,
    """
    f_fft = fftpack.fftshift(fftpack.fftn(f))
    g_fft = fftpack.fftshift(fftpack.fftn(g))
    return fftpack.fftshift(fftpack.ifftn(fftpack.ifftshift(f_fft/g_fft))) 

def _slp_filter(phase, cutoff, rows, cols, x_size, y_size, params):
    """
    Function to perform spatial low pass filter
    """
    cx = np.floor(cols/2)
    cy = np.floor(rows/2)
    # fft for the input image
    imf = fftshift(fft2(phase))
    # calculate distance
    distfact = 1.0e3  # to convert into meters
    [xx, yy] = np.meshgrid(range(cols), range(rows))
    xx = (xx - cx) * x_size  # these are in meters as x_size in meters
    yy = (yy - cy) * y_size
    dist = np.sqrt(xx ** 2 + yy ** 2)/distfact  # km

    if params[cf.SLPF_METHOD] == 1:  # butterworth low pass filter
        H = 1. / (1 + ((dist / cutoff) ** (2 * params[cf.SLPF_ORDER])))
    else:  # Gaussian low pass filter
        H = np.exp(-(dist ** 2) / (2 * cutoff ** 2))
    outf = imf * H
    out = np.real(ifft2(ifftshift(outf)))
    out[np.isnan(phase)] = np.nan
    return out  # out is units of phase, i.e. mm


# TODO: use tiles here and distribute amongst processes 

def _get_autogrid(phase):
    """
    Helper function to assist with memory re-allocation during FFT calculation
    """
    autocorr_grid = _calc_autoc_grid(phase)
    nzc = np.sum(np.sum(phase != 0))
    autocorr_grid = fftshift(real(autocorr_grid)) / nzc
    return autocorr_grid 

def two_point_correlation_fft(im):
    r"""
    Calculates the two-point correlation function using fourier transforms

    Parameters
    ----------
    im : ND-array
        The image of the void space on which the 2-point correlation is desired

    Returns
    -------
    result : named_tuple
        A tuple containing the x and y data for plotting the two-point
        correlation function, using the *args feature of matplotlib's plot
        function.  The x array is the distances between points and the y array
        is corresponding probabilities that points of a given distance both
        lie in the void space.

    Notes
    -----
    The fourier transform approach utilizes the fact that the autocorrelation
    function is the inverse FT of the power spectrum density.
    For background read the Scipy fftpack docs and for a good explanation see:
    http://www.ucl.ac.uk/~ucapikr/projects/KamilaSuankulova_BSc_Project.pdf
    """
    # Calculate half lengths of the image
    hls = (np.ceil(np.shape(im))/2).astype(int)
    # Fourier Transform and shift image
    F = sp_ft.ifftshift(sp_ft.fftn(sp_ft.fftshift(im)))
    # Compute Power Spectrum
    P = np.absolute(F**2)
    # Auto-correlation is inverse of Power Spectrum
    autoc = np.absolute(sp_ft.ifftshift(sp_ft.ifftn(sp_ft.fftshift(P))))
    tpcf = _radial_profile(autoc, r_max=np.min(hls))
    return tpcf 

def mfcc(s,fs, nfiltbank):
  
    #divide into segments of 25 ms with overlap of 10ms
    nSamples = np.int32(0.025*fs)
    overlap = np.int32(0.01*fs)
    nFrames = np.int32(np.ceil(len(s)/(nSamples-overlap)))
    #zero padding to make signal length long enough to have nFrames
    padding = ((nSamples-overlap)*nFrames) - len(s)
    if padding > 0:
        signal = np.append(s, np.zeros(padding))
    else:
        signal = s
    segment = np.empty((nSamples, nFrames))
    start = 0
    for i in range(nFrames):
        segment[:,i] = signal[start:start+nSamples]
        start = (nSamples-overlap)*i
    
    #compute periodogram
    nfft = 512
    periodogram = np.empty((nFrames, int(nfft/2 + 1)))
    for i in range(nFrames):
        x = segment[:,i] * hamming(nSamples)
        spectrum = fftshift(fft(x,nfft))
        periodogram[i,:] = abs(spectrum[int(nfft/2-1):])/nSamples
        
    #calculating mfccs    
    fbank = mel_filterbank(nfft, nfiltbank, fs)
    #nfiltbank MFCCs for each frame
    mel_coeff = np.empty((nfiltbank,nFrames))
    for i in range(nfiltbank):
        for k in range(nFrames):
            mel_coeff[i,k] = np.sum(periodogram[k,:]*fbank[:,i])
            
    mel_coeff = np.log10(mel_coeff)
    mel_coeff = dct(mel_coeff)
    #exclude 0th order coefficient (much larger than others)
    mel_coeff[0,:]= np.zeros(nFrames)
    return mel_coeff 

def correlation_2D(image):
    """
    #TODO document normalization output in units

    :param image: 2d image
    :return: 2d fourier transform
    """
    # Take the fourier transform of the image.
    F1 = fftpack.fft2(image)

    # Now shift the quadrants around so that low spatial frequencies are in
    # the center of the 2D fourier transformed image.
    F2 = fftpack.fftshift(F1)
    return np.abs(F2) 

def SO3_FFT_synthesize(f_hat):
    """
    Perform the inverse (spectral to spatial) SO(3) Fourier transform.

    :param f_hat: a list of matrices of with shapes [1x1, 3x3, 5x5, ..., 2 L_max + 1 x 2 L_max + 1]
    """
    F = wigner_d_transform_synthesis(f_hat)

    # The rest of the SO(3) FFT is just a standard torus FFT
    F = fftshift(F, axes=(0, 2))
    f = ifft2(F, axes=(0, 2))

    b = len(f_hat)
    return f * (2 * b) ** 2 

def generate_fractal_surface(self, G):
        """Generate a 2D array with a fractal distribution.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        if self.xs == self.xf:
            surfacedims = (self.ny, self.nz)
        elif self.ys == self.yf:
            surfacedims = (self.nx, self.nz)
        elif self.zs == self.zf:
            surfacedims = (self.nx, self.ny)

        self.fractalsurface = np.zeros(surfacedims, dtype=complextype)

        # Positional vector at centre of array, scaled by weighting
        v1 = np.array([self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1] * (surfacedims[1]) / 2])

        # 2D array of random numbers to be convolved with the fractal function
        R = np.random.RandomState(self.seed)
        A = R.randn(surfacedims[0], surfacedims[1])

        # 2D FFT
        A = fftpack.fftn(A)
        # Shift the zero frequency component to the centre of the array
        A = fftpack.fftshift(A)

        # Generate fractal
        generate_fractal2D(surfacedims[0], surfacedims[1], G.nthreads, self.b, self.weighting, v1, A, self.fractalsurface)

        # Shift the zero frequency component to start of the array
        self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
        # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT
        self.fractalsurface = np.real(fftpack.ifftn(self.fractalsurface))
        # Scale the fractal volume according to requested range
        fractalmin = np.amin(self.fractalsurface)
        fractalmax = np.amax(self.fractalsurface)
        fractalrange = fractalmax - fractalmin
        self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) \
            + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin 

def test_inverse(self):
        for n in [1,4,9,100,211]:
            x = random.random((n,))
            assert_array_almost_equal(ifftshift(fftshift(x)),x) 

def test_definition(self):
        x = [0,1,2,3,4,-4,-3,-2,-1]
        y = [-4,-3,-2,-1,0,1,2,3,4]
        assert_array_almost_equal(fftshift(x),y)
        assert_array_almost_equal(ifftshift(y),x)
        x = [0,1,2,3,4,-5,-4,-3,-2,-1]
        y = [-5,-4,-3,-2,-1,0,1,2,3,4]
        assert_array_almost_equal(fftshift(x),y)
        assert_array_almost_equal(ifftshift(y),x) 

def test_inverse(self):
        for n in [1,4,9,100,211]:
            x = random((n,))
            assert_array_almost_equal(ifftshift(fftshift(x)),x) 

def test_definition(self):
        x = [0,1,2,3,4,-4,-3,-2,-1]
        y = [-4,-3,-2,-1,0,1,2,3,4]
        assert_array_almost_equal(fftshift(x),y)
        assert_array_almost_equal(ifftshift(y),x)
        x = [0,1,2,3,4,-5,-4,-3,-2,-1]
        y = [-5,-4,-3,-2,-1,0,1,2,3,4]
        assert_array_almost_equal(fftshift(x),y)
        assert_array_almost_equal(ifftshift(y),x) 

def convolve(f, g):
    """
    FFT based convolution

    :param f: array
    :param g: array
    :return: array, (f * g)[n]
    """
    f_fft = fftpack.fftshift(fftpack.fftn(f))
    g_fft = fftpack.fftshift(fftpack.fftn(g))
    return fftpack.fftshift(fftpack.ifftn(fftpack.ifftshift(f_fft*g_fft))) 

def fourier_fft(fEM, time, freq, ftarg):
    r"""Fourier Transform using the Fast Fourier Transform.

    The function is called from one of the modelling routines in
    :mod:`empymod.model`. Consult these modelling routines for a description of
    the input and output parameters.

    Returns
    -------
    tEM : array
        Returns time-domain EM response of `fEM` for given `time`.

    conv : bool
        Only relevant for QWE/QUAD.

    """
    # Get ftarg values
    dfreq = ftarg['dfreq']
    nfreq = ftarg['nfreq']
    ntot = ftarg['ntot']
    pts_per_dec = ftarg['pts_per_dec']

    # If pts_per_dec, we have first to interpolate fEM to required freqs
    if pts_per_dec:
        sfEMr = iuSpline(np.log(freq), fEM.real)
        sfEMi = iuSpline(np.log(freq), fEM.imag)
        freq = np.arange(1, nfreq+1)*dfreq
        fEM = sfEMr(np.log(freq)) + 1j*sfEMi(np.log(freq))

    # Pad the frequency result
    fEM = np.pad(fEM, (0, ntot-nfreq), 'linear_ramp')

    # Carry out FFT
    ifftEM = fftpack.ifft(np.r_[fEM[1:], 0, fEM[::-1].conj()]).real
    stEM = 2*ntot*fftpack.fftshift(ifftEM*dfreq, 0)

    # Interpolate in time domain
    dt = 1/(2*ntot*dfreq)
    ifEM = iuSpline(np.linspace(-ntot, ntot-1, 2*ntot)*dt, stEM)
    tEM = ifEM(time)/2*np.pi  # (Multiplication of 2/pi in model.tem)

    # Return the electromagnetic time domain field
    # (Second argument is only for QWE)
    return tEM, True


# 3. Utilities 

def stft(signal, fs, nfft, overlap):

    #plotting time domain signal
    plt.figure(1)
    t = np.arange(0,len(signal)/fs, 1/fs)
    plt.plot(t,signal)
    plt.axis(xmax = 1)
    plt.xlabel('Time in seconds')
    plt.ylabel('Amplitude')
    plt.title('Speech signal')
    
    if not np.log2(nfft).is_integer():
        nfft = nearestPow2(nfft)
    slength = len(signal)
    hop_size = np.int32(overlap * nfft)
    nFrames = int(np.round(len(signal)/(nfft-hop_size)))
    #zero padding to make signal length long enough to have nFrames
    signal = np.append(signal, np.zeros(nfft))
    STFT = np.empty((nfft, nFrames))
    segment = np.zeros(nfft)
    start = 0
    
    for n in range(nFrames):
        segment = signal[start:start+nfft] * hann(nfft) 
        padded_seg = np.append(segment,np.zeros(nfft))
        spec = fftshift(fft(padded_seg))
        spec = spec[len(spec)/2:]
        spec = abs(spec)/max(abs(spec))
        powerspec = 20*np.log10(spec)
        STFT[:,n] = powerspec
        start = start + nfft - hop_size 
        
    #plot spectrogram
    plt.figure(2)
    freq = (fs/(2*nfft)) * np.arange(0,nfft,1)
    time = np.arange(0,nFrames)*(slength/(fs*nFrames))
    plt.imshow(STFT, extent = [0,max(time),0,max(freq)],origin='lower', cmap='jet', interpolation='nearest', aspect='auto')
    plt.ylabel('Frequency in Hz')
    plt.xlabel('Time in seconds')
    plt.axis([0,max(time),0,np.max(freq)])
    plt.title('Spectrogram of speech')
    plt.show()
    return (STFT, time, freq) 

def SO3_FFT_analyze(f):
    """
    Compute the complex SO(3) Fourier transform of f.

    The standard way to define the FT is:
    \hat{f}^l_mn = (2 J + 1)/(8 pi^2) *
       int_0^2pi da int_0^pi db sin(b) int_0^2pi dc f(a,b,c) D^{l*}_mn(a,b,c)

    The normalizing constant comes about because:
    int_SO(3) D^*(g) D(g) dg = 8 pi^2 / (2 J + 1)
    where D is any Wigner D function D^l_mn. Note that the factor 8 pi^2 (the volume of SO(3))
    goes away if we integrate with the normalized Haar measure.

    This function computes the FT using the normalized D functions:
    \tilde{D} = 1/2pi sqrt((2J+1)/2) D
    where D are the rotation matrices in the basis of complex, seismology-normalized, centered spherical harmonics.

    Hence, this function computes:
    \hat{f}^l_mn = \int_SO(3) f(g) \tilde{D}^{l*}_mn(g) dg

    So that the FT of f = \tilde{D}^l_mn is 1 at (l,m,n) (and zero elsewhere).

    Args:
      f: an array of shape (2B, 2B, 2B), where B is the bandwidth.

    Returns:
      f_hat: the Fourier transform of f. A list of length B,
      where entry l contains an 2l+1 by 2l+1 array containing the projections
      of f onto matrix elements of the l-th irreducible representation of SO(3).

    Main source:
    SOFT: SO(3) Fourier Transforms
    Peter J. Kostelec and Daniel N. Rockmore

    Further information:
    Generalized FFTs-a survey of some recent results
    Maslen & Rockmore

    Engineering Applications of Noncommutative Harmonic Analysis.
    9.5 - Sampling and FFT for SO(3) and SU(2)
    G.S. Chrikjian, A.B. Kyatkin
    """
    assert f.shape[0] == f.shape[1]
    assert f.shape[1] == f.shape[2]
    assert f.shape[0] % 2 == 0

    # First, FFT along the alpha and gamma axes (axis 0 and 2, respectively)
    F = fft2(f, axes=(0, 2))
    F = fftshift(F, axes=(0, 2))

    # Then, perform the Wigner-d transform
    return wigner_d_transform_analysis(F) 

def generate_fractal_volume(self, G):
        """Generate a 3D volume with a fractal distribution.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        # Scale filter according to size of fractal volume
        if self.nx == 1:
            filterscaling = np.amin(np.array([self.ny, self.nz])) / np.array([self.ny, self.nz])
            filterscaling = np.insert(filterscaling, 0, 1)
        elif self.ny == 1:
            filterscaling = np.amin(np.array([self.nx, self.nz])) / np.array([self.nx, self.nz])
            filterscaling = np.insert(filterscaling, 1, 1)
        elif self.nz == 1:
            filterscaling = np.amin(np.array([self.nx, self.ny])) / np.array([self.nx, self.ny])
            filterscaling = np.insert(filterscaling, 2, 1)
        else:
            filterscaling = np.amin(np.array([self.nx, self.ny, self.nz])) / np.array([self.nx, self.ny, self.nz])

        # Adjust weighting to account for filter scaling
        self.weighting = np.multiply(self.weighting, filterscaling)

        self.fractalvolume = np.zeros((self.nx, self.ny, self.nz), dtype=complextype)

        # Positional vector at centre of array, scaled by weighting
        v1 = np.array([self.weighting[0] * self.nx / 2, self.weighting[1] * self.ny / 2, self.weighting[2] * self.nz / 2])

        # 3D array of random numbers to be convolved with the fractal function
        R = np.random.RandomState(self.seed)
        A = R.randn(self.nx, self.ny, self.nz)

        # 3D FFT
        A = fftpack.fftn(A)
        # Shift the zero frequency component to the centre of the array
        A = fftpack.fftshift(A)

        # Generate fractal
        generate_fractal3D(self.nx, self.ny, self.nz, G.nthreads, self.b, self.weighting, v1, A, self.fractalvolume)

        # Shift the zero frequency component to the start of the array
        self.fractalvolume = fftpack.ifftshift(self.fractalvolume)
        # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT
        self.fractalvolume = np.real(fftpack.ifftn(self.fractalvolume))
        # Bin fractal values
        bins = np.linspace(np.amin(self.fractalvolume), np.amax(self.fractalvolume), self.nbins)
        for j in range(self.ny):
            for k in range(self.nz):
                self.fractalvolume[:, j, k] = np.digitize(self.fractalvolume[:, j, k], bins, right=True)