def kappa_to_gamma(kappa,dt1,dt2=None):
    """
    simple application of Kaiser-Squires (1995) kernel in fourier
    space to convert complex shear to complex convergence: imaginary
    part of convergence is B-mode.
    """
    if not dt2:
        dt2 = dt1
    N1,N2 = kappa.shape

    #convert angles from arcminutes to radians
    dt1 = dt1 * np.pi / 180. / 60.
    dt2 = dt2 * np.pi / 180. / 60.

    #compute k values corresponding to field size
    dk1 = np.pi / N1 / dt1
    dk2 = np.pi / N2 / dt2

    k1 = fftpack.ifftshift( dk1 * (np.arange(2*N1)-N1) )
    k2 = fftpack.ifftshift( dk2 * (np.arange(2*N2)-N2) )

    ipart,rpart = np.meshgrid(k2,k1)
    k = rpart + 1j*ipart

    #compute Kaiser-Squires kernel on this grid. Eq. 43 p. 329
    fourier_map = np.conj( KS_kernel(k) )

    #compute Fourier transform of the kappa
    kappa_fft = fftpack.fft2( kappa, (2*N1,2*N2) )

    gamma_fft = kappa_fft * fourier_map

    gamma = fftpack.ifft2(gamma_fft)[:N1,:N2]

    return gamma

def spec2grid(sfield):
    """
    Transform one frame of SQG model
    output to a grided (physical) representation.  Assumes 'sfield'
    to be up-half plane, and specifies lower half plane by conjugate
    sym (since physical field is assumed real-valued).  Input field
    should have dimensions  (...,kmax+1,2*kmax+1,nz), where
    kmax=2^n-1, hence physical resolution will be 2^(n+1) x 2^(n+1).
    NOTE: top row of the input field corresponds to ky = 0, the
    kx<0 part is NOT assumed a priori to be conjugate- symmetric
    with the kx>0 part.  NOTE: grid2spec(spec2grid(fk)) = fk.
    OPTIONAL: da = true pads input with 0s before transfoming to
    gridspace, for dealiased products.  Default is da = false.
    
    Args:
        sfield: complex spectrum field with shape (t(optional), ky, kx, z(optional))
    """        
    if not _is_single_layer(sfield):
        hres = sfield.shape[-2] + 1
        fk = fullspec(sfield)
        fk = fftpack.ifftshift(fk, axes=(-2,-3))
        return hres*hres*np.real(fftpack.ifft2(fk, axes=(-2,-3)))
    else:
        hres = sfield.shape[-1] + 1
        fk = fullspec(sfield, True)
        fk = fftpack.ifftshift(fk, axes=(-1,-2))
        return hres*hres*np.real(fftpack.ifft2(fk, axes=(-1,-2)))

def frankotchellappa(dzdx, dzdy):
    rows, cols = dzdx.shape

    # The following sets up matrices specifying frequencies in the x and y
    # directions corresponding to the Fourier transforms of the gradient
    # data.  They range from -0.5 cycles/pixel to + 0.5 cycles/pixel.
    # The scaling of this is irrelevant as long as it represents a full
    # circle domain. This is functionally equivalent to any constant * pi
    pi_over_2 = np.pi / 2.0
    row_grid = np.linspace(-pi_over_2, pi_over_2, rows)
    col_grid = np.linspace(-pi_over_2, pi_over_2, cols)
    wy, wx = np.meshgrid(row_grid, col_grid, indexing='ij')

    # Quadrant shift to put zero frequency at the appropriate edge
    wx = ifftshift(wx)
    wy = ifftshift(wy)

    # Fourier transforms of gradients
    DZDX = fft2(dzdx)
    DZDY = fft2(dzdy)

    # Integrate in the frequency domain by phase shifting by pi/2 and
    # weighting the Fourier coefficients by their frequencies in x and y and
    # then dividing by the squared frequency
    denom = (wx ** 2 + wy ** 2)
    Z = (-1j * wx * DZDX - 1j * wy * DZDY) / denom
    Z = np.nan_to_num(Z)

    return np.real(ifft2(Z))

def convolve2D_analytic(theta,F):
    func = KS_kernel_fourier
    #func = gaussian_fourier

    assert theta.shape == F.shape
    
    dtheta1 = theta[1,1].real - theta[0,0].real
    dtheta2 = theta[1,1].imag - theta[0,0].imag
    
    N1,N2 = theta.shape

    dell1 = numpy.pi / N1 / dtheta1
    dell2 = numpy.pi / N2 / dtheta2
    
    ell1 = fftpack.ifftshift( dell1 * (numpy.arange(2*N1)-N1) )
    ell2 = fftpack.ifftshift( dell2 * (numpy.arange(2*N2)-N2) )
    
    ell = numpy.zeros((2*N1,2*N2),dtype=complex)
    ell += ell1.reshape((2*N1,1))
    ell += 1j * ell2

    F_fft = fftpack.fft2(F,(2*N1,2*N2) )
    F_fft *= func( ell )

    out = fftpack.ifft2(F_fft)

    return out[:N1,:N2]

def convolve2D_fft(theta,F):
    func = KS_kernel_real
    #func = gaussian_real

    assert theta.shape == F.shape

    N1,N2 = theta.shape
    
    dtheta1 = theta[1,1].real - theta[0,0].real
    theta1 = dtheta1*(numpy.arange(2*N1)-N1)
    theta1 = fftpack.ifftshift(theta1)
    
    dtheta2 = theta[1,1].imag - theta[0,0].imag
    theta2 = dtheta2*(numpy.arange(2*N2)-N2)
    theta2 = fftpack.ifftshift(theta2)
    
    theta_kernel = numpy.zeros((2*N1,2*N2),dtype=complex)
    theta_kernel += theta1.reshape((2*N1,1))
    theta_kernel += 1j*theta2

    kernel = func(theta_kernel)
    
    dA = dtheta1 * dtheta2

    F_fft = fftpack.fft2(F, (2*N1,2*N2) ) * dA
    F_fft *= fftpack.fft2(kernel,(2*N1,2*N2) ) 
    
    out = fftpack.ifft2(F_fft)
    
    return out[:N1,:N2]

def gamma_to_kappa(shear,dt1,dt2=None):
    """
    simple application of Kaiser-Squires (1995) kernel in fourier
    space to convert complex shear to complex convergence: imaginary
    part of convergence is B-mode.
    """
    if not dt2:
        dt2 = dt1
    N1,N2 = shear.shape

    #convert angles from arcminutes to radians
    dt1 = dt1 * numpy.pi / 180. / 60.
    dt2 = dt2 * numpy.pi / 180. / 60.

    #compute k values corresponding to field size
    dk1 = numpy.pi / N1 / dt1
    dk2 = numpy.pi / N2 / dt2

    k1 = fftpack.ifftshift( dk1 * (numpy.arange(2*N1)-N1) )
    k2 = fftpack.ifftshift( dk2 * (numpy.arange(2*N2)-N2) )

    ipart,rpart = numpy.meshgrid(k2,k1)
    k = rpart + 1j*ipart

    #compute (inverse) Kaiser-Squires kernel on this grid
    fourier_map = numpy.conj( KS_kernel(-k) )
    
    #compute Fourier transform of the shear
    gamma_fft = fftpack.fft2( shear, (2*N1,2*N2) )

    kappa_fft = fourier_map * gamma_fft

    kappa = fftpack.ifft2(kappa_fft)[:N1,:N2]

    return kappa

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 fftPropagate(field, grid, propDistance):
    '''Propagates a sampled 1D field along the optical axis.
    
    fftPropagate propagates a sampled 1D field a distance L by computing the
    field's angular spectrum, multiplying each spectral component by the
    propagation kernel exp(j * 2 * pi * fx * L / wavelength), and then
    recomibining the propagated spectral components. The angular spectrum is
    computed using a FFT.
    
    Parameters
    ----------
    field        : 1D array of complex
        The sampled field to propagate.
    grid         : Grid
        The grid on which the sampled field lies.
    propDistance : float
        The distance to propagate the field in the same physical units as the
        grid.
    
    '''
    scalingFactor = (grid.physicalSize / (grid.gridSize - 1))
    F             = scalingFactor * fftshift(fft(ifftshift(field)))
    
    # Compute the z-component of the wavevector
    # Adding 0j ensures that numpy.sqrt returns complex numbers
    kz = 2 * np.pi * np.sqrt(1 - (grid.pfX * grid.wavelength)**2 + 0j) / grid.wavelength
    
    # Propagate the field's spectral components
    Fprop = F * np.exp(1j * kz * propDistance)
    
    # Recombine the spectral components
    fieldProp = fftshift(ifft(ifftshift(Fprop))) / scalingFactor
    
    return fieldProp

 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 preWhitenCube(**kwargs):
	'''
	Pre-whitenening using noise estimates from a cube taken from the difference map.
	Returns a the pre-whitened volume and various spectra. (Alp Kucukelbir, 2013)

	'''
	print '\n= Pre-whitening the Cubes'
	tStart = time()

	n             = kwargs.get('n', 0)
	vxSize        = kwargs.get('vxSize', 0)
	elbowAngstrom = kwargs.get('elbowAngstrom', 0)
	rampWeight    = kwargs.get('rampWeight',1.0)
	dataF         = kwargs.get('dataF', 0)
	dataBGF       = kwargs.get('dataBGF', 0)
	dataBGSpect   = kwargs.get('dataBGSpect', 0)

	epsilon = 1e-10

	pWfilter = createPreWhiteningFilter(n             = n,
										spectrum      = dataBGSpect,
										elbowAngstrom = elbowAngstrom,
										rampWeight    = rampWeight,
										vxSize        = vxSize)

	# Apply the pre-whitening filter to the inside cube
	dataF       = np.multiply(pWfilter['pWfilter'],dataF)

	dataPWFabs  = np.abs(dataF)
	dataPWFabs  = dataPWFabs-np.min(dataPWFabs)
	dataPWFabs  = dataPWFabs/np.max(dataPWFabs)
	dataPWSpect = sphericalAverage(dataPWFabs**2) + epsilon

	dataPW = np.real(fftpack.ifftn(fftpack.ifftshift(dataF)))
	del dataF

	# Apply the pre-whitening filter to the outside cube
	dataBGF       = np.multiply(pWfilter['pWfilter'],dataBGF)

	dataPWBGFabs  = np.abs(dataBGF)
	dataPWBGFabs  = dataPWBGFabs-np.min(dataPWBGFabs)
	dataPWBGFabs  = dataPWBGFabs/np.max(dataPWBGFabs)
	dataPWBGSpect = sphericalAverage(dataPWBGFabs**2) + epsilon

	dataBGPW = np.real(fftpack.ifftn(fftpack.ifftshift(dataBGF)))
	del dataBGF

	m, s = divmod(time() - tStart, 60)
	print "  :: Time elapsed: %d minutes and %.2f seconds" % (m, s)

	return {'dataPW':dataPW, 'dataBGPW':dataBGPW, 'dataPWSpect': dataPWSpect, 'dataPWBGSpect': dataPWBGSpect, 'peval': pWfilter['peval'], 'pcoef': pWfilter['pcoef'] }

def filtergrid(rows, cols):

    # Set up u1 and u2 matrices with ranges normalised to +/- 0.5
    u1, u2 = np.meshgrid(np.linspace(-0.5, 0.5, cols, endpoint=(cols % 2)),
                         np.linspace(-0.5, 0.5, rows, endpoint=(rows % 2)),
                         sparse=True)

    # Quadrant shift to put 0 frequency at the top left corner
    u1 = ifftshift(u1)
    u2 = ifftshift(u2)

    # Compute frequency values as a radius from centre (but quadrant shifted)
    radius = np.sqrt(u1 * u1 + u2 * u2)

    return radius, u1, u2

    def Afunc(self, f):
        fs = reshape(f, (self.height, self.width, self.depth))

        F = fftn(fs)

        d_1 = ifftshift(ifftn(F*self.H));
        
        d_e = ifftshift(ifftn(F*self.He));
        
        d_e2 = ifftshift(ifftn(F*self.He2));
        
        d = (1.5*real(d_1) + 2*real(d_e*self.e1) + 0.5*real(d_e2*self.e2))
        
        d = real(d);
        return ravel(d)

    def Ahfunc(self, f):
        fs = reshape(f, (self.height, self.width, self.depth))

        F = fftn(fs)

        d_1 = ifftshift(ifftn(F*self.Ht));

        d_e = ifftshift(ifftn(F*self.Het));

        d_e2 = ifftshift(ifftn(F*self.He2t));

        d = (1.5*d_1 + 2*real(d_e*exp(1j*self.alpha)) + 0.5*real(d_e2*exp(2*1j*self.alpha)));
        
        d = real(d);
        return ravel(d)

def plot_pulse(A_t,A_w,t,w, l0 = 1550 * nm, t_zoom = ps, l_zoom = 10*nm):
        ## Fix maximum
        if 'A_max' not in plot_pulse.__dict__:
            plot_pulse.A_max = amax(A_t)
            plot_pulse.A_w_max = amax(A_w)

        w0 = 2 * pi * C_SPEED / l0

        ## Plot Time domain
        fig = pl.gcf()
        fig.add_subplot(211)
        pl.plot(t / t_zoom, absolute(A_t/plot_pulse.A_max), hold = False)
        pl.axis([amin(t) / t_zoom*0.4,amax(t) / t_zoom*0.4, 0, 1.1])
        pl.xlabel(r'$time\ (%s s)$'%units[t_zoom])

        atten_win = 0.01
        npoints = len(w)
        apod = arange(npoints)
        apod = exp(-apod/(npoints*atten_win)) + exp(-(npoints-apod)/(npoints*atten_win))
        apod = ifftshift(apod)

        ## Plot Freq domain
        fig.add_subplot(212)
        pl.plot((2 * pi * C_SPEED)/(w+w0)/nm, log10(absolute(A_w)**2/absolute(plot_pulse.A_w_max)**2) * 10, hold=False)
        #pl.plot((2 * pi * C_SPEED)/(w+w0) / Q_('um'), log10(absolute(apod)**2/absolute(amax(apod))**2 ) *10, hold=True)
        #pl.semilogy()
        pl.axis([(l0-l_zoom)/Q_('um'), (l0+l_zoom)/nm, -60, 5])
        pl.xlabel(r'$wavelength\ (\mu m)$')
        pl.ylabel(r'$spectrum (db)$')

        pl.show()

 def invertPsd1d(self, psdx=None, psd1d=None, phasespec=None, seed=None):
     """Convert a 1d PSD, generate a phase spectrum (or use user-supplied values) into a 2d PSD (psd2dI)."""
     # Converting the 2d PSD into a 1d PSD is a lossy process, and then there is additional randomness
     #  added when the phase spectrum is not the same as the original phase spectrum, so this may or may not
     #  look that much like the original image (unless you keep the phases, doesn't look like image). 
     # 'Swing' the 1d PSD across the whole fov.         
     if psd1d == None:
         psd1d = self.psd1d
         xr = self.rfreq / self.xfreqscale
     else:
         if psdx == None:
             xr = numpy.arange(0, len(psd1d), 1.0)
     # Resample into even bins (definitely necessarily if using psd1d).
     xrange = numpy.arange(0, numpy.sqrt(self.xcen**2 + self.ycen**2)+1.0, 1.0)
     psd1d = numpy.interp(xrange, xr, psd1d, right=psd1d[len(psd1d)-1])
     # Calculate radii - distance from center.
     rad = numpy.hypot((self.yy-self.ycen), (self.xx-self.xcen))
     # Calculate the PSD2D from the 1d value.
     self.psd2dI = numpy.interp(rad.flatten(), xrange, psd1d)
     self.psd2dI = self.psd2dI.reshape(self.ny, self.nx)
     if phasespec == None:
         self._makeRandomPhases(seed=seed)
     else:
         self.phasespecI = phasespec
     if not(self.shift):
         # The 1d PSD is centered, so the 2d PSD will be 'shifted' here. 
         self.psd2dI = fftpack.ifftshift(self.psd2dI)
     return

 def invertPsd2d(self, useI=False, usePhasespec=True, seed=None):
     """Convert the 2d PSD and phase spec into an FFT image (FftI). """
     # The PHASEs of the FFT are encoded in the phasespec ('where things are')
     # The AMPLITUDE of the FFT is encoded in the PSD ('how bright things are' .. also radial scale)
     # amp = sqrt(| R(uv)^2 + I(u,v)^2|) == length of 'z' (z = x + iy, in fourier image)
     # phase = arctan(y/x)
     if useI:
         psd2d = self.psd2dI
         phasespec = self.phasespecI
     else:
         psd2d = self.psd2d
         phasespec = self.phasespec
     if self.shift:
         amp = numpy.sqrt(fftpack.ifftshift(psd2d))
     else:
         amp = numpy.sqrt(psd2d)
     # Can override using own phasespec for random (if coming in here to use 2d PSD for image)
     if not(usePhasespec):
         self._makeRandomPhases(seed=seed)
         phasespec = self.phasespecI
     # Shift doesn't matter for phases, because 'unshifted' it above, before calculating phase.
     x = numpy.cos(phasespec) * amp
     y = numpy.sin(phasespec) * amp
     self.fimageI = x + 1j*y
     if self.shift:
         self.fimageI = fftpack.fftshift(self.fimageI)
     return

	def _FFT_raw_fired(self):
		fft_shift = fft.fft2(Data.TrA_Data)
		Data.FFT = fft.ifftshift(fft_shift)
		plt.figure()
		plt.contourf(np.log(np.abs(Data.FFT)**2),cmap=plt.cm.Greys_r)
		plt.title('FFT of raw')
		plt.show()

 def pdf(self):
     """ Applies the 3D FFT in the q-space grid to generate
     the DSI diffusion propagator, remove the background noise with a
     hard threshold and then deconvolve the propagator with the 
     Lucy-Richardson deconvolution algorithm
     """
     values = self.data
     #create the signal volume
     Sq = np.zeros((self.qgrid_sz, self.qgrid_sz, self.qgrid_sz))
     #fill q-space
     for i in range(self.dn):
         qx, qy, qz = self.model.qgrid[i]
         Sq[qx, qy, qz] += values[i]
     #get deconvolution PSF
     DSID_PSF = self.model.cache_get('deconv_psf', key=self.model.gtab)
     if DSID_PSF is None:
         DSID_PSF = gen_PSF(self.model.qgrid, self.qgrid_sz, 
                            self.qgrid_sz, self.qgrid_sz)
     self.model.cache_set('deconv_psf', self.model.gtab, DSID_PSF)
     #apply fourier transform
     Pr = fftshift(np.abs(np.real(fftn(ifftshift(Sq), 
                   3 * (self.qgrid_sz, )))))
     #threshold propagator
     Pr = threshold_propagator(Pr)
     #apply LR deconvolution
     Pr = LR_deconv(Pr, DSID_PSF, 5, 2)
     return Pr

 def initialize_time_freq(self):
     if self.thread_initialize_tfr is not None or self.is_computing.any():
         # needd to come back later ...
         if not self.timer_back_initialize.isActive():
             self.timer_back_initialize.start()
         return
     
     # create self.params_time_freq
     p = self.params_time_freq = { }
     for param in self.paramTimeFreq.children():
         self.params_time_freq[param.name()] = param.value()
     
     
     # we take sampling_rate = f_stop*4 or (original sampling_rate)
     if p['f_stop']*4 < self.global_sampling_rate:
         p['sampling_rate'] = p['f_stop']*4
     else:
         p['sampling_rate']  = self.global_sampling_rate
     self.factor = p['sampling_rate']/self.global_sampling_rate # this compensate unddersampling in FFT.
     
     self.xsize2 = self.xsize
     self.len_wavelet = int(self.xsize2*p['sampling_rate'])
     self.win = fftpack.ifftshift(np.hamming(self.len_wavelet))
     self.thread_initialize_tfr = ThreadInitializeWavelet(len_wavelet = self.len_wavelet, 
                                                         params_time_freq = p, parent = self )
     self.thread_initialize_tfr.finished.connect(self.initialize_tfr_done)
     self.thread_initialize_tfr.start()

def spectralWhitening(data, sr=None, smoothi=None, freq_domain=False, apply_filter=None):
    """
    Apply spectral whitening to data.
    
    sr: sampling rate (only needed for smoothing)
    smoothi: None or int
    Data is divided by its smoothed (Default: None) amplitude spectrum.
    """
    if freq_domain:
        mask = False
        spec = data
    else:
        mask = np.ma.getmask(data)
        N = len(data)
        nfft = nextpow2(N)
        spec = fft(data, nfft)
    #df = sr/N
    spec_ampl = np.sqrt(np.abs(np.multiply(spec, np.conjugate(spec))))
    if isinstance(smoothi, basestring) and isnumber(smoothi) and smoothi > 0:
        smoothi = int(smoothi * N / sr)
        spec /= ifftshift(smooth(fftshift(spec_ampl), smoothi))
    else:
        spec /= spec_ampl
    if apply_filter is not None:
        spec *= filterResp(*apply_filter, sr=sr, N=len(spec), whole=True)[1]
    if freq_domain:
        return spec
    else:
        ret = np.real(ifft(spec, nfft)[:N])
        if USE_FFTW3:
            ret = ret.copy()
        return fillArray(ret, mask=mask, fill_value=0.)