Python numpy.fft.fftfreq() 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 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 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 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 test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def spd(self, npos):
        '''raw spectral density, returns Fourier transform

        n is number of points in positive spectrum, the actual number of points
        is twice as large. different from other spd methods with fft
        '''
        n = npos
        w = fft.fftfreq(2*n) * 2 * np.pi
        hw = self.fftarma(2*n)  #not sure, need to check normalization
        #return (hw*hw.conj()).real[n//2-1:]  * 0.5 / np.pi #doesn't show in plot
        return (hw*hw.conj()).real * 0.5 / np.pi, w 

def spdshift(self, n):
        '''power spectral density using fftshift

        currently returns two-sided according to fft frequencies, use first half
        '''
        #size = s1+s2-1
        mapadded = self.padarr(self.ma, n)
        arpadded = self.padarr(self.ar, n)
        hw = fft.fft(fft.fftshift(mapadded)) / fft.fft(fft.fftshift(arpadded))
        #return np.abs(spd)[n//2-1:]
        w = fft.fftfreq(n) * 2 * np.pi
        wslice = slice(n//2-1, None, None)
        #return (hw*hw.conj()).real[wslice], w[wslice]
        return (hw*hw.conj()).real, w 

def spddirect(self, n):
        '''power spectral density using padding to length n done by fft

        currently returns two-sided according to fft frequencies, use first half
        '''
        #size = s1+s2-1
        #abs looks wrong
        hw = fft.fft(self.ma, n) / fft.fft(self.ar, n)
        w = fft.fftfreq(n) * 2 * np.pi
        wslice = slice(None, n//2, None)
        #return (np.abs(hw)**2)[wslice], w[wslice]
        return (np.abs(hw)**2) * 0.5/np.pi, w 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def spd(self, npos):
        '''raw spectral density, returns Fourier transform

        n is number of points in positive spectrum, the actual number of points
        is twice as large. different from other spd methods with fft
        '''
        n = npos
        w = fft.fftfreq(2*n) * 2 * np.pi
        hw = self.fftarma(2*n)  #not sure, need to check normalization
        #return (hw*hw.conj()).real[n//2-1:]  * 0.5 / np.pi #doesn't show in plot
        return (hw*hw.conj()).real * 0.5 / np.pi, w 

def spdshift(self, n):
        '''power spectral density using fftshift

        currently returns two-sided according to fft frequencies, use first half
        '''
        #size = s1+s2-1
        mapadded = self.padarr(self.ma, n)
        arpadded = self.padarr(self.ar, n)
        hw = fft.fft(fft.fftshift(mapadded)) / fft.fft(fft.fftshift(arpadded))
        #return np.abs(spd)[n//2-1:]
        w = fft.fftfreq(n) * 2 * np.pi
        wslice = slice(n//2-1, None, None)
        #return (hw*hw.conj()).real[wslice], w[wslice]
        return (hw*hw.conj()).real, w 

def spddirect(self, n):
        '''power spectral density using padding to length n done by fft

        currently returns two-sided according to fft frequencies, use first half
        '''
        #size = s1+s2-1
        #abs looks wrong
        hw = fft.fft(self.ma, n) / fft.fft(self.ar, n)
        w = fft.fftfreq(n) * 2 * np.pi
        wslice = slice(None, n//2, None)
        #return (np.abs(hw)**2)[wslice], w[wslice]
        return (np.abs(hw)**2) * 0.5/np.pi, w 

def frequencies(self):
        return fftfreq(self.n_fft_samples, 1.0 / self.sampling_frequency) 

def main(args):
    pyfftw.interfaces.cache.enable()

    refmap = mrc.read(args.key, compat="relion")
    df = star.parse_star(args.input, keep_index=False)
    star.augment_star_ucsf(df)
    refmap_ft = vop.vol_ft(refmap, threads=args.threads)

    apix = star.calculate_apix(df)
    sz = refmap_ft.shape[0] // 2 - 1
    sx, sy = np.meshgrid(rfftfreq(sz), fftfreq(sz))
    s = np.sqrt(sx ** 2 + sy ** 2)
    r = s * sz
    r = np.round(r).astype(np.int64)
    r[r > sz // 2] = sz // 2 + 1
    a = np.arctan2(sy, sx)

    def1 = df["rlnDefocusU"].values
    def2 = df["rlnDefocusV"].values
    angast = df["rlnDefocusAngle"].values
    phase = df["rlnPhaseShift"].values
    kv = df["rlnVoltage"].values
    ac = df["rlnAmplitudeContrast"].values
    cs = df["rlnSphericalAberration"].values
    xshift = df["rlnOriginX"].values
    yshift = df["rlnOriginY"].values

    score = np.zeros(df.shape[0])

    # TODO parallelize
    for i, row in df.iterrows():
        xcor = particle_xcorr(row, refmap_ft)

    if args.top is None:
        args.top = df.shape[0]

    top = df.iloc[np.argsort(score)][:args.top]
    star.simplify_star_ucsf(top)
    star.write_star(args.output, top)
    return 0 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def test_definition(self):
        x = [0, 1, 2, 3, 4, -4, -3, -2, -1]
        assert_array_almost_equal(9*fft.fftfreq(9), x)
        assert_array_almost_equal(9*pi*fft.fftfreq(9, pi), x)
        x = [0, 1, 2, 3, 4, -5, -4, -3, -2, -1]
        assert_array_almost_equal(10*fft.fftfreq(10), x)
        assert_array_almost_equal(10*pi*fft.fftfreq(10, pi), x) 

def component_viewer(output, tr=2.0):
    ''' This a function to interactively view the results of a decomposition analysis

    Args:
        output: (dict) output dictionary from running Brain_data.decompose()
        tr: (float) repetition time of data
    '''

    if ipywidgets is None:
        raise ImportError(
            "ipywidgets is required for interactive plotting. Please install this package manually or install nltools with optional arguments: pip install 'nltools[interactive_plots]'"
        )

    def component_inspector(component, threshold):
        '''This a function to be used with ipywidgets to interactively view a decomposition analysis

            Make sure you have tr and output assigned to variables.

            Example:

                from ipywidgets import BoundedFloatText, BoundedIntText
                from ipywidgets import interact

                tr = 2.4
                output = data_filtered_smoothed.decompose(algorithm='ica', n_components=30, axis='images', whiten=True)

                interact(component_inspector, component=BoundedIntText(description='Component', value=0, min=0, max=len(output['components'])-1),
                      threshold=BoundedFloatText(description='Threshold', value=2.0, min=0, max=4, step=.1))

        '''
        _, ax = plt.subplots(nrows=3, figsize=(12,8))
        thresholded = (output['components'][component] - output['components'][component].mean())*(1/output['components'][component].std())
        thresholded.data[np.abs(thresholded.data) <= threshold] = 0
        plot_stat_map(thresholded.to_nifti(), cut_coords=range(-40, 70, 10),
                      display_mode='z', black_bg=True, colorbar=True, annotate=False,
                      draw_cross=False, axes=ax[0])
        if isinstance(output['decomposition_object'], (sklearn.decomposition.PCA)):
            var_exp = output['decomposition_object'].explained_variance_ratio_[component]
            ax[0].set_title(f"Component: {component}/{len(output['components'])}, Variance Explained: {var_exp:2.2}", fontsize=18)
        else:
            ax[0].set_title(f"Component: {component}/{len(output['components'])}", fontsize=18)

        ax[1].plot(output['weights'][:, component], linewidth=2, color='red')
        ax[1].set_ylabel('Intensity (AU)', fontsize=18)
        ax[1].set_title(f'Timecourse (TR={tr})', fontsize=16)
        y = fft(output['weights'][:, component])
        f = fftfreq(len(y), d=tr)
        ax[2].plot(f[f > 0], np.abs(y)[f > 0]**2)
        ax[2].set_ylabel('Power', fontsize=18)
        ax[2].set_xlabel('Frequency (Hz)', fontsize=16)

    ipywidgets.interact(component_inspector, component=ipywidgets.BoundedIntText(description='Component', value=0, min=0, max=len(output['components'])-1),
              threshold=ipywidgets.BoundedFloatText(description='Threshold', value=2.0, min=0, max=4, step=.1))