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

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 find_foci(evt, type,key,type2,key2,minPhase=-500000, maxPhase=500000, steps=101, field_of_view_rad=100, wavelength=1.053, CCD_S_DIST=0.375, PX_SIZE=75e-6):
    img = evt[type][key].data
    centroids = evt[type2][key2].data

    Nfoci = centroids.shape[0]
    Xrange, Yrange = img.shape
    Npixel = field_of_view_rad
    
    p = numpy.linspace(-Xrange/2, Xrange/2-1, Xrange)
    q = numpy.linspace(-Yrange/2, Yrange/2-1, Yrange)
    pp, qq = numpy.meshgrid(p, q)
   
    phase_matrix = (2*numpy.pi/wavelength)*numpy.sqrt(1-((PX_SIZE/CCD_S_DIST)**2)*(qq**2 + pp**2))
    prop_length = numpy.linspace(minPhase, maxPhase, steps)
    
    variance = numpy.zeros([steps, Nfoci])
    # shift stuff for performance reasons
    img_shifted = fftshift(img)
    phase_matrix_shifted = fftshift(phase_matrix)
    
    for idx, phase in enumerate(prop_length):
        
        img_propagated = img_shifted * numpy.exp(1.j*phase*phase_matrix_shifted)
        recon = fftshift(ifft2(img_propagated))
        
        for CC in numpy.arange(Nfoci):
            centerx, centery = centroids[CC, :]
            ###print  centerx, centery
            reconcut = numpy.abs(recon[numpy.max([0, centerx-Npixel-1]).astype(int): numpy.min([Xrange-1, centerx+Npixel]).astype(int), numpy.max([0, centery-Npixel-1]).astype(int): numpy.min([Yrange-1, centery+Npixel]).astype(int)])
            variance[idx, CC] = reconcut.var()
    
    focus_distance = numpy.zeros(Nfoci)
    CC_size = numpy.zeros(Nfoci)
    focused_CC = numpy.zeros(4*Npixel**2 * Nfoci).reshape(Nfoci, 2*Npixel, 2*Npixel)
    
    for CC in numpy.arange(Nfoci):
        ind_max = numpy.argmax(variance[:, CC])
        tmp = variance[:, CC]
        # get max which is not at border
        loc_max_bool = numpy.r_[True, tmp[1:] > tmp[:-1]] & numpy.r_[tmp[:-1] > tmp[1:], True]
        loc_max_bool[0] = False
        loc_max_bool[-1] = False
        ind_max = numpy.argmax(tmp*loc_max_bool)
        
        focus_distance[CC] = prop_length[ind_max]
        img_propagated = img_shifted * numpy.exp(1.j * focus_distance[CC]  * phase_matrix_shifted)
        recon = fftshift(ifft2(img_propagated))
        
        centerx, centery = centroids[CC, :]
        
        reconcut = numpy.real(recon[numpy.max([0, centerx-Npixel]).astype(int): numpy.min([Xrange-1, centerx+Npixel]).astype(int), numpy.max([0, centery-Npixel]).astype(int): numpy.min([Yrange-1, centery+Npixel]).astype(int)])
        focused_CC[CC, 0:reconcut.shape[0], 0:reconcut.shape[1]] = reconcut
        CC_size[CC] = numpy.sum(get_CC_size(reconcut))
    
    if len(focused_CC):
        add_record(evt["analysis"], "analysis", "focused_CC", focused_CC[0])
        add_record(evt["analysis"], "analysis", "focus distance", focus_distance)
        add_record(evt["analysis"], "analysis", "CC_size", CC_size)
        add_record(evt["analysis"], "analysis", "propagation length", prop_length)

    def fourierTransform(self, fromPos, toPos, only = []):
        self.checkToPos(toPos)
        if len(only) > 0:
            self.allFid[toPos] = np.array([fftshift(fft(self.allFid[fromPos][fidIndex])) for fidIndex in only])
        else:
            self.allFid[toPos] = np.array([fftshift(fft(fid)) for fid in self.allFid[fromPos]])

        self.frequency = np.linspace(-self.sweepWidthTD2/2,self.sweepWidthTD2/2,len(self.allFid[fromPos][0]))

def plot_fft_shifted(x, y, *args, **kwargs):
    import pylab
    symmetric = False
    if "symmetric" in kwargs:
        symmetric = kwargs.pop("symmetric")
    if symmetric and np.mod(x.shape[0],2)==0:
        pylab.plot(fftp.fftshift(x)[1:], fftp.fftshift(y)[1:], *args, **kwargs)
    else:
        pylab.plot(fftp.fftshift(x), fftp.fftshift(y), *args, **kwargs)

def f2(file1, file2):
  from scipy.fftpack import fftshift
  (imgarr1,w,h) = pgm.pgmread(file1)
  (imgarr2,w,h) = pgm.pgmread(file2)
  imgarr1 = np.float32(imgarr1)
  imgarr2 = np.float32(imgarr2)
  print imgarr1.shape, imgarr2.shape
  #
  dftArr1 = dft_2d_from_1d(imgarr1)
  dftArr2 = dft_2d_from_1d(imgarr2)
  #
  absArr1 = np.abs(dftArr1)
  absArr2 = np.abs(dftArr2)
  phaseArr1 = np.angle(dftArr1)
  phaseArr2 = np.angle(dftArr2)
  #
  newArr1 = absArr1 * np.exp(1j * phaseArr2)
  newArr2 = absArr2 * np.exp(1j * phaseArr1)
  #print 'Error : ', np.sum((newArr1 - dftArr1)**2)
  #print 'Error : ', np.sum((newArr2 - dftArr2)**2)
  #
  idftArr1 = idft_2d_from_1d(newArr1)
  idftArr2 = idft_2d_from_1d(newArr2)
  idftArr1 = np.round(idftArr1.real)
  idftArr2 = np.round(idftArr2.real)
  #
  dftArrAbs1 = np.abs(fftshift(dftArr1))
  dftArrAbs2 = np.abs(fftshift(dftArr2))
  #
  plt.figure()
  plt.subplot(121)
  plt.title(file1)
  plt.imshow(imgarr1, cmap=cm.gray)
  #plt.colorbar()
  plt.subplot(122)
  plt.title(file2)
  plt.imshow(imgarr2, cmap=cm.gray)
  #plt.colorbar()
  #
  plt.figure()
  plt.subplot(121)
  plt.title('DFT of ' + file1)
  plt.imshow(np.int32(np.log(dftArrAbs1)), cmap=cm.gray)
  plt.colorbar()
  plt.subplot(122)
  plt.title('DFT of ' + file2)
  plt.imshow(np.int32(np.log(dftArrAbs2)), cmap=cm.gray)
  plt.colorbar()
  #
  plt.figure()
  plt.subplot(121)
  plt.title('IDFT1')
  plt.imshow(np.int32(idftArr1), cmap=cm.gray)
  #plt.colorbar()
  plt.subplot(122)
  plt.title('IDFT2')
  plt.imshow(np.int32(idftArr2), cmap=cm.gray)