Python numpy.fft.fft2() Examples

The following are 23 code examples for showing how to use numpy.fft.fft2(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy.fft , or try the search function .

Example 1
Project: opticspy   Author: Sterncat   File: zernike.py    License: MIT License 6 votes vote down vote up
def ptf(self):
		"""
		Phase transfer function
		"""
		PSF = self.__psfcaculator__()
		PTF = __fftshift__(__fft2__(PSF))
		PTF = __np__.angle(PTF)
		b = 400
		R = (200)**2
		for i in range(b):
			for j in range(b):
				if (i-b/2)**2+(j-b/2)**2>R:
					PTF[i][j] = 0
		__plt__.imshow(abs(PTF),cmap=__cm__.rainbow)
		__plt__.colorbar()
		__plt__.show()
		return 0 
Example 2
Project: prnu-python   Author: polimi-ispl   File: functions.py    License: MIT License 6 votes vote down vote up
def wiener_dft(im: np.ndarray, sigma: float) -> np.ndarray:
    """
    Adaptive Wiener filter applied to the 2D FFT of the image
    :param im: multidimensional array
    :param sigma: estimated noise power
    :return: filtered version of input im
    """
    noise_var = sigma ** 2
    h, w = im.shape

    im_noise_fft = fft2(im)
    im_noise_fft_mag = np.abs(im_noise_fft / (h * w) ** .5)

    im_noise_fft_mag_noise = wiener_adaptive(im_noise_fft_mag, noise_var)

    zeros_y, zeros_x = np.nonzero(im_noise_fft_mag == 0)

    im_noise_fft_mag[zeros_y, zeros_x] = 1
    im_noise_fft_mag_noise[zeros_y, zeros_x] = 0

    im_noise_fft_filt = im_noise_fft * im_noise_fft_mag_noise / im_noise_fft_mag
    im_noise_filt = np.real(ifft2(im_noise_fft_filt))

    return im_noise_filt.astype(np.float32) 
Example 3
Project: prnu-python   Author: polimi-ispl   File: functions.py    License: MIT License 6 votes vote down vote up
def crosscorr_2d(k1: np.ndarray, k2: np.ndarray) -> np.ndarray:
    """
    PRNU 2D cross-correlation
    :param k1: 2D matrix of size (h1,w1)
    :param k2: 2D matrix of size (h2,w2)
    :return: 2D matrix of size (max(h1,h2),max(w1,w2))
    """
    assert (k1.ndim == 2)
    assert (k2.ndim == 2)

    max_height = max(k1.shape[0], k2.shape[0])
    max_width = max(k1.shape[1], k2.shape[1])

    k1 -= k1.flatten().mean()
    k2 -= k2.flatten().mean()

    k1 = np.pad(k1, [(0, max_height - k1.shape[0]), (0, max_width - k1.shape[1])], mode='constant', constant_values=0)
    k2 = np.pad(k2, [(0, max_height - k2.shape[0]), (0, max_width - k2.shape[1])], mode='constant', constant_values=0)

    k1_fft = fft2(k1, )
    k2_fft = fft2(np.rot90(k2, 2), )

    return np.real(ifft2(k1_fft * k2_fft)).astype(np.float32) 
Example 4
Project: pytorch_fft   Author: locuslab   File: test.py    License: Apache License 2.0 6 votes vote down vote up
def test_expand(): 
    X = torch.randn(2,2,4,4).cuda().double()
    zeros = torch.zeros(2,2,4,4).cuda().double()
    r1, r2 = cfft.rfft2(X)
    c1, c2 = cfft.fft2(X, zeros)
    assert np.allclose(cfft.expand(r1).cpu().numpy(), c1.cpu().numpy())
    assert np.allclose(cfft.expand(r2, imag=True).cpu().numpy(), c2.cpu().numpy())
    r1, r2 = cfft.rfft3(X)
    c1, c2 = cfft.fft3(X, zeros)
    assert np.allclose(cfft.expand(r1).cpu().numpy(), c1.cpu().numpy())
    assert np.allclose(cfft.expand(r2, imag=True).cpu().numpy(), c2.cpu().numpy())

    X = torch.randn(2,2,5,5).cuda().double()
    zeros = torch.zeros(2,2,5,5).cuda().double()
    r1, r2 = cfft.rfft3(X)
    c1, c2 = cfft.fft3(X, zeros)
    assert np.allclose(cfft.expand(r1, odd=True).cpu().numpy(), c1.cpu().numpy())
    assert np.allclose(cfft.expand(r2, imag=True, odd=True).cpu().numpy(), c2.cpu().numpy()) 
Example 5
Project: pysteps   Author: pySTEPS   File: fft.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
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 
Example 6
Project: pysteps   Author: pySTEPS   File: fft.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
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 
Example 7
Project: lie_learn   Author: AMLab-Amsterdam   File: T2FFT.py    License: MIT License 5 votes vote down vote up
def analyze(f, axes=(0, 1)):
        """
        Compute the Fourier Transform of the discretely sampled function f : T^2 -> C.

        Let f : T^2 -> C be a band-limited function on the torus.
        The samples f(theta_k, phi_l) correspond to points on a regular grid on the circle,
        as returned by spaces.T1.linspace:
        theta_k = phi_k = 2 pi k / N
        for k = 0, ..., N - 1 and l = 0, ..., N - 1

        This function computes
        \hat{f}_n = (1/N) \sum_{k=0}^{N-1} f(theta_k) e^{-i n theta_k}
        which, if f has band-limit less than N, is equal to:
        \hat{f}_n = \int_0^{2pi} f(theta) e^{-i n theta} dtheta / 2pi,
                  = <f(theta), e^{i n theta}>
        where dtheta / 2pi is the normalized Haar measure on T^1, and < , > denotes the inner product on Hilbert space,
        with respect to which this transform is unitary.

        The range of frequencies n is -floor(N/2) <= n <= ceil(N/2) - 1

        :param f:
        :param axis:
        :return:
        """
        # The numpy FFT returns coefficients in a different order than we want them,
        # and using a different normalization.
        f_hat = fft2(f, axes=axes)
        f_hat = fftshift(f_hat, axes=axes)
        size = np.prod([f.shape[ax] for ax in axes])
        return f_hat / size 
Example 8
Project: opticspy   Author: Sterncat   File: zernike_rec.py    License: MIT License 5 votes vote down vote up
def __psfcaculator__(self,lambda_1=632*10**(-9),z=0.1):
		"""
		height: Exit pupil height
		width: Exit pupil width
		z: Distance from exit pupil to image plane
		"""
		a = self.__a__
		b = __sqrt__(1-a**2)
		l1 = 100;
		x1 = __np__.linspace(-a, a, l1)
		y1 = __np__.linspace(-b, b, l1)
		[X,Y] = __np__.meshgrid(x1,y1)
		Z = __zernikecartesian__(self.__coefficients__,a,X,Y)
		d = 400 # background
		A = __np__.zeros([d,d])
		A[d//2-l1//2+1:d//2+l1//2+1,d//2-l1//2+1:d//2+l1//2+1] = Z
		# fig = __plt__.figure()
		# __plt__.imshow(A)
		# __plt__.colorbar()
		# __plt__.show()
		abbe = __np__.exp(-1j*2*__np__.pi*A)
		for i in range(len(abbe)):
			for j in range(len(abbe)):
				if abbe[i][j]==1:
					abbe[i][j]=0
		PSF = __fftshift__(__fft2__(__fftshift__(abbe)))**2
		PSF = PSF/PSF.max()
		return PSF 
Example 9
Project: opticspy   Author: Sterncat   File: zernike_rec.py    License: MIT License 5 votes vote down vote up
def mtf(self,lambda_1=632*10**(-9),z=0.1,matrix = False):
		"""
		Modulate Transfer function
		"""
		PSF = self.__psfcaculator__(lambda_1=lambda_1,z=z)
		MTF = __fftshift__(__fft2__(PSF))
		MTF = MTF/MTF.max()
		fig = __plt__.figure(figsize=(9, 6), dpi=80)
		__plt__.imshow(abs(MTF),cmap=__cm__.bwr)
		__plt__.colorbar()
		__plt__.show()
		if matrix == True:
			return MTF
		else:
			return 0 
Example 10
Project: opticspy   Author: Sterncat   File: zernike_rec.py    License: MIT License 5 votes vote down vote up
def ptf(self):
		"""
		Phase transfer function
		"""
		PSF = self.__psfcaculator__()
		PTF = __fftshift__(__fft2__(PSF))
		PTF = __np__.angle(PTF)
		l1 = 100
		d = 400
		A = __np__.zeros([d,d])
		A[d//2-l1//2+1:d//2+l1//2+1,d//2-l1//2+1:d//2+l1//2+1] = PTF[d//2-l1//2+1:d//2+l1//2+1,d//2-l1//2+1:d//2+l1//2+1]
		__plt__.imshow(abs(A),cmap=__cm__.rainbow)
		__plt__.colorbar()
		__plt__.show()
		return 0 
Example 11
Project: opticspy   Author: Sterncat   File: zernike.py    License: MIT License 5 votes vote down vote up
def __psfcaculator__(self,r=1,lambda_1=632*10**(-9),z=0.1):
		"""
		pupil: Exit pupil diameter
		z: Distance from exit pupil to image plane
		r: pupil radius, in unit of lambda
		"""
		pupil = l1 = 200 # exit pupil sample points
		x = __np__.linspace(-r, r, l1)
		[X,Y] = __np__.meshgrid(x,x)
		Z = __interferometer__.__zernikecartesian__(self.__coefficients__,X,Y)
		for i in range(len(Z)):
			for j in range(len(Z)):
				if x[i]**2+x[j]**2>r**2:
					Z[i][j] = 0
		d = 400 # background
		A = __np__.zeros([d,d])
		A[d//2-l1//2+1:d//2+l1//2+1,d//2-l1//2+1:d//2+l1//2+1] = Z
		axis_1 = d//pupil*r
		fig = __plt__.figure()
		# ax = fig.gca()
		# __plt__.imshow(A,extent=[-axis_1,axis_1,-axis_1,axis_1],cmap=__cm__.RdYlGn)
		# ax.set_xlabel('mm',fontsize=14)
		# __plt__.colorbar()
		# __plt__.show()

		abbe = __np__.exp(-1j*2*__np__.pi*A)
		for i in range(len(abbe)):
			for j in range(len(abbe)):
				if abbe[i][j]==1:
					abbe[i][j]=0
		PSF = __fftshift__(__fft2__(__fftshift__(abbe)))**2
		PSF = PSF/PSF.max()
		return PSF 
Example 12
Project: opticspy   Author: Sterncat   File: zernike.py    License: MIT License 5 votes vote down vote up
def otf(self,r=1,lambda_1=632*10**(-9),z=0.1):
		PSF = self.__psfcaculator__(r=r,lambda_1=lambda_1,z=z)
		OTF = __fftshift__(__fft2__(PSF))
		return 0 
Example 13
Project: ProxImaL   Author: comp-imaging   File: utils.py    License: MIT License 5 votes vote down vote up
def fftd(I, dims=None):

    # Compute fft
    if dims is None:
        X = fftn(I)
    elif dims == 2:
        X = fft2(I, axes=(0, 1))
    else:
        X = fftn(I, axes=tuple(range(dims)))

    return X 
Example 14
Project: ray-optics   Author: mjhoptics   File: analyses.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def calc_psf(wavefront, ndim, maxdim):
    """Calculate the point spread function of wavefront W.

    Args:
        wavefront: ndim x ndim Numpy array of wavefront errors. No data
                   condition is indicated by nan
        ndim: The sampling across the wavefront
        maxdim: The total width of the sampling grid

    Returns: AP, the PSF of the input wavefront
    """
    maxdim_by_2 = maxdim//2
    W = np.zeros([maxdim, maxdim])
    nd2 = ndim//2
    W[maxdim_by_2-(nd2-1):maxdim_by_2+(nd2+1),
      maxdim_by_2-(nd2-1):maxdim_by_2+(nd2+1)] = np.nan_to_num(wavefront)

    phase = np.exp(1j*2*np.pi*W)

    for i in range(len(phase)):
        for j in range(len(phase)):
            if phase[i][j] == 1:
                phase[i][j] = 0

    AP = abs(fftshift(fft2(fftshift(phase))))**2
    AP_max = np.nanmax(AP)
    AP = AP/AP_max
    return AP 
Example 15
Project: picasso   Author: jungmannlab   File: imageprocess.py    License: MIT License 5 votes vote down vote up
def xcorr(imageA, imageB):
    FimageA = _fft.fft2(imageA)
    CFimageB = _np.conj(_fft.fft2(imageB))
    return _fft.fftshift(
        _np.real(_fft.ifft2((FimageA * CFimageB)))
    ) / _np.sqrt(imageA.size) 
Example 16
Project: tierpsy-tracker   Author: ver228   File: helper.py    License: MIT License 5 votes vote down vote up
def fft_convolve2d(x,y):
    """ 2D convolution, using FFT"""
    fr = fft2(x)
    fr2 = fft2(y)
    cc = np.real(ifft2(fr*fr2))
    cc = fftshift(cc)
    return cc 
Example 17
Project: game-of-life   Author: thearn   File: lib.py    License: Apache License 2.0 5 votes vote down vote up
def fft_convolve2d(x,y):
    """
    2D convolution, using FFT
    """
    fr = fft2(x)
    fr2 = fft2(np.flipud(np.fliplr(y)))
    m,n = fr.shape
    cc = np.real(ifft2(fr*fr2))
    cc = np.roll(cc, - int(m / 2) + 1, axis=0)
    cc = np.roll(cc, - int(n / 2) + 1, axis=1)
    return cc 
Example 18
Project: odl   Author: odlgroup   File: shearlab_operator.py    License: Mozilla Public License 2.0 5 votes vote down vote up
def sheardec2D(X, shearletsystem):
    """Shearlet Decomposition function."""
    coeffs = np.zeros(shearletsystem.shearlets.shape, dtype=complex)
    Xfreq = fftshift(fft2(ifftshift(X)))
    for i in range(shearletsystem.nShearlets):
        coeffs[:, :, i] = fftshift(ifft2(ifftshift(Xfreq * np.conj(
                                   shearletsystem.shearlets[:, :, i]))))
    return coeffs.real 
Example 19
Project: odl   Author: odlgroup   File: shearlab_operator.py    License: Mozilla Public License 2.0 5 votes vote down vote up
def shearrec2D(coeffs, shearletsystem):
    """Shearlet Recovery function."""
    X = np.zeros(coeffs.shape[:2], dtype=complex)
    for i in range(shearletsystem.nShearlets):
        X = X + fftshift(fft2(
            ifftshift(coeffs[:, :, i]))) * shearletsystem.shearlets[:, :, i]
    return (fftshift(ifft2(ifftshift((
            X / shearletsystem.dualFrameWeights))))).real 
Example 20
Project: odl   Author: odlgroup   File: shearlab_operator.py    License: Mozilla Public License 2.0 5 votes vote down vote up
def sheardecadjoint2D(coeffs, shearletsystem):
    """Shearlet Decomposition adjoint function."""
    X = np.zeros(coeffs.shape[:2], dtype=complex)
    for i in range(shearletsystem.nShearlets):
        X = X + fftshift(fft2(
            ifftshift(coeffs[:, :, i]))) * np.conj(
            shearletsystem.shearlets[:, :, i])
    return (fftshift(ifft2(ifftshift(
            X / shearletsystem.dualFrameWeights)))).real 
Example 21
Project: odl   Author: odlgroup   File: shearlab_operator.py    License: Mozilla Public License 2.0 5 votes vote down vote up
def shearrecadjoint2D(X, shearletsystem):
    """Shearlet Recovery adjoint function."""
    coeffs = np.zeros(shearletsystem.shearlets.shape, dtype=complex)
    Xfreq = fftshift(fft2(ifftshift(X)))
    for i in range(shearletsystem.nShearlets):
        coeffs[:, :, i] = fftshift(ifft2(ifftshift(
            Xfreq * shearletsystem.shearlets[:, :, i])))
    return coeffs.real 
Example 22
Project: pyoptools   Author: cihologramas   File: gs.py    License: GNU General Public License v3.0 4 votes vote down vote up
def gs(idata,itera=10, ia=None):
    """Gerchberg-Saxton algorithm to calculate DOEs
    
    Calculates the phase distribution in a object plane to obtain an 
    specific amplitude distribution in the target plane. It uses a 
    FFT to calculate the field propagation.
    The wavefront at the DOE plane is assumed as a plane wave.
    
    **ARGUMENTS:**
    
        ========== ========================================================
        idata      numpy array containing the target amplitude distribution 
        itera      Maximum number of iterations
        ia         Illumination amplitude at the hologram plane if not given
                   it is assumed to be a constant amplitude with a value
                   of 1. If given it should be an array with the same shape
                   of idata
        ========== ========================================================
    """
    
    if ia==None:
        inpa=ones(idata.shape)
    else:
        inpa=ia
    
    assert idata.shape==inpa.shape, "ia and idata must have the same dimensions"
    
    fdata=fftshift(fft2(ifftshift(idata)))
    e=1000
    ea=1000
    
    for i in range (itera):
        fdata=exp(1.j*angle(fdata))*inpa
        
        rdata=ifftshift(ifft2(fftshift(fdata)))
        e= (abs(rdata)-idata).std()
        if e>ea: 
            break
        ea=e
        rdata=exp(1.j*angle(rdata))*(idata)
        fdata=fftshift(fft2(ifftshift(rdata)))        
    
    fdata=exp(1.j*angle(fdata))
    return fdata*inpa 
Example 23
Project: pyoptools   Author: cihologramas   File: gs.py    License: GNU General Public License v3.0 4 votes vote down vote up
def gs_mod(idata,itera=10,osize=256):
    """Modiffied Gerchberg-Saxton algorithm to calculate DOEs
    
    Calculates the phase distribution in a object plane to obtain an 
    specific amplitude distribution in the target plane. It uses a 
    FFT to calculate the field propagation.
    The wavefront at the DOE plane is assumed as a plane wave.
    This algorithm leaves a window around the image plane to allow the 
    noise to move there. It only optimises the center of the image.
    
    **ARGUMENTS:**
    
        ========== ======================================================
        idata      numpy array containing the target amplitude distribution 
        itera      Maximum number of iterations
        osize      Size of the center of the image to be optimized
                   It should be smaller than the image itself.
        ========== ======================================================
    """
    M,N=idata.shape
    cut=osize//2
    
    
    zone=zeros_like(idata)
    zone[M/2-cut:M/2+cut,N/2-cut:N/2+cut]=1
    zone=zone.astype(bool)

    mask=exp(2.j*pi*random(idata.shape))
    mask[zone]=0
    
    #~ imshow(abs(mask)),colorbar()
    
    fdata=fftshift(fft2(ifftshift(idata+mask))) #Nota, colocar esta mascara es muy importante, por que si no  no converge tan rapido
    
    e=1000
    ea=1000
    for i in range (itera):
        fdata=exp(1.j*angle(fdata))

        rdata=ifftshift(ifft2(fftshift(fdata)))
        #~ e= (abs(rdata[zone])-idata[zone]).std()
        #~ if e>ea: 
           #~ 
            #~ break
        ea=e
        rdata[zone]=exp(1.j*angle(rdata[zone]))*(idata[zone])        
        fdata=fftshift(fft2(ifftshift(rdata)))   
    fdata=exp(1.j*angle(fdata))
    return fdata