def _UpsampledDFT(data, nor, noc, precision=1, roff=0, coff=0):
    """
    Upsampled DFT by matrix multiplies.
    data (numpy.array): 2d array
    nor, noc (ints): Number of pixels in the output upsampled DFT, in units
    of upsampled pixels
    precision (int): Calculate drift within 1/precision of a pixel
    roff, coff (ints): Row and column offsets, allow to shift the output array
                    to a region of interest on the DFT
    returns (tuple of floats): Drift in pixels
    """
    z = 1j  # imaginary unit
    nr, nc = data.shape

    # Compute kernels and obtain DFT by matrix products
    kernc = numpy.power(math.e, (-z * 2 * math.pi / (nc * precision)) *
                                ((fft.ifftshift(arange(0, nc))[:, None]).T - nc // 2) *
                                (arange(0, noc) - coff)[:, None]
                       )

    kernr = numpy.power(math.e, (-z * 2 * math.pi / (nr * precision)) *
                                (fft.ifftshift(arange(0, nr))[:, None] - nr // 2) *
                                ((arange(0, nor)[:, None]).T - roff)
                       )

    return numpy.dot(numpy.dot((kernr.transpose()), data), kernc.transpose())

def compare_interpolated_spectrum():
    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111)

    out = fft(ifftshift(f_full))
    freqs = fftfreq(len(f_full), d=0.01) # spacing, Ang
    sfreqs = fftshift(freqs)
    taper = gauss_taper(freqs, sigma=0.0496) #Ang, corresponds to 2.89 km/s at 5150A.
    tout = out * taper

    ax.plot(sfreqs, fftshift(tout))

    wl_h, fl_h = np.abs(np.load("PH6.8kms_0.01ang.npy"))
    wl_l, fl_l = np.abs(np.load("PH2.5kms.npy"))

    #end edges
    wl_he = wl_h[200:-200]
    fl_he = fl_h[200:-200]
    interp = Sinc_w(wl_l, fl_l, a=5, window='kaiser')
    fl_hi = interp(wl_he)

    d = wl_he[1] - wl_he[0]
    out = fft(ifftshift(fl_hi))
    freqs = fftfreq(len(out), d=d)
    ax.plot(fftshift(freqs), fftshift(out))

    plt.show()

def whiten_olsh_lee_inner(image, f_0=None, central_clip=(None, None), normalize_pre=True, normalize_post=True,
                          no_filter=False):
    height, width = image.shape
    assert height % 2 == 0 and width % 2 == 0, "image must have even size!"
    if normalize_pre:
        image = image - image.mean()  # I personally think this is useless, since rho will make (0,0) freq compoenent 0.
        std_im = image.std(ddof=1)
        assert std_im != 0, "constant image unsupported!"
        image /= std_im

    fx, fy = np.meshgrid(np.arange(-height / 2, height / 2), np.arange(-width / 2, width / 2), indexing='ij')
    rho = np.sqrt(fx * fx + fy * fy)
    if f_0 is None:
        f_0 = 0.4 * (height + width) / 2
    filt = rho * np.exp(-((rho / f_0) ** 4))

    im_f = fft2(image)
    if not no_filter:
        fft_filtered_old = im_f * ifftshift(filt)
    else:  # hack to only lower frequency response.
        print('no real filtering!')
        fft_filtered_old = im_f
    fft_filtered_old = fftshift(fft_filtered_old)
    if central_clip != (None, None):
        fft_filtered_old = fft_filtered_old[height // 2 - central_clip[0] // 2:height // 2 + central_clip[0] // 2,
                           width // 2 - central_clip[1] // 2:width // 2 + central_clip[1] // 2]
    im_out = np.real(ifft2(ifftshift(fft_filtered_old)))
    # I believe since the rho at the (0,0) frequency part is zero, then the whole image should be zero as well.
    # so explicit DC removing is useless.
    if normalize_post:
        assert abs(im_out.mean()) < 1e-6  # should be extremely small.
        std_im_out = im_out.std(ddof=1)
    else:
        std_im_out = 1
    return im_out / std_im_out

    def _simulate_image(self):
        """
        Generates the fake output.
        """
        with self._acquisition_init_lock:
            pos = self.align.position.value
            logging.debug("Simulating image shift by %s", pos)
            ac, bc = pos.get("a"), pos.get("b")
            ang = math.radians(135)
            # AB->XY
            xc = -(ac * math.sin(ang) + bc * math.cos(ang))
            yc = -(ac * math.cos(ang) - bc * math.sin(ang))
            pixelSize = self.fake_img.metadata[model.MD_PIXEL_SIZE]
            self.fake_img.metadata[model.MD_ACQ_DATE] = time.time()
            x_pxs = xc / pixelSize[0]
            y_pxs = yc / pixelSize[1]

            # Image shifted based on LensAligner position
            z = 1j  # imaginary unit
            self.deltar = x_pxs
            self.deltac = y_pxs
            nr, nc = self.fake_img.shape
            array_nr = numpy.arange(-numpy.fix(nr / 2), numpy.ceil(nr / 2))
            array_nc = numpy.arange(-numpy.fix(nc / 2), numpy.ceil(nc / 2))
            Nr = fft.ifftshift(array_nr)
            Nc = fft.ifftshift(array_nc)
            [Nc, Nr] = numpy.meshgrid(Nc, Nr)
            sim_img = fft.ifft2(fft.fft2(self.fake_img) * numpy.power(math.e,
                            z * 2 * math.pi * (self.deltar * Nr / nr + self.deltac * Nc / nc)))
            output = model.DataArray(abs(sim_img), self.fake_img.metadata)
            return output

def odddftups(inp,nor=None,noc=None,usfac=1,roff=0,coff=0):
    from numpy.fft import ifftshift
    from numpy import pi,newaxis,floor

    nr,nc=np.shape(inp);

    # Set defaults
    if noc is None: noc=nc;
    if nor is None: nor=nr;

    if nr % 2 == 1:
        oddr = True
        nrnew = nr+1
    else:
        oddr = False
    if nr % 2 == 1:
        oddr = True
        nrnew = nr+1
    else:
        oddr = False

    # Compute kernels and obtain DFT by matrix products
    term1c = ( ifftshift(np.arange(nc) - floor(nc/2)).T[:,newaxis] )
    term2c = ( np.arange(noc) - coff  )[newaxis,:]
    kernc=np.exp((-1j*2*pi/(nc*usfac))*term1c*term2c);
    term1r = ( np.arange(nor).T - roff )[:,newaxis]
    term2r = ( ifftshift(np.arange(nr)) - floor(nr/2) )[newaxis,:]
    kernr=np.exp((-1j*2*pi/(nr*usfac))*term1r*term2r);
    #kernc=exp((-i*2*pi/(nc*usfac))*( ifftshift([0:nc-1]).' - floor(nc/2) )*( [0:noc-1] - coff ));
    #kernr=exp((-i*2*pi/(nr*usfac))*( [0:nor-1].' - roff )*( ifftshift([0:nr-1]) - floor(nr/2)  ));
    out=np.dot(np.dot(kernr,inp),kernc);
    #return np.roll(np.roll(out,+1,axis=0),+1,axis=1)
    return out 

 def test_axes_keyword(self):
     freqs = [[ 0,  1,  2], [ 3,  4, -4], [-3, -2, -1]]
     shifted = [[-1, -3, -2], [ 2,  0,  1], [-4,  3,  4]]
     assert_array_almost_equal(fftshift(freqs, axes=(0, 1)), shifted)
     assert_array_almost_equal(fftshift(freqs, axes=0), fftshift(freqs, axes=(0,)))
     assert_array_almost_equal(ifftshift(shifted, axes=(0, 1)), freqs)
     assert_array_almost_equal(ifftshift(shifted, axes=0), ifftshift(shifted, axes=(0,)))

    def fft_g2r(self, fg, fg_ishifted=False):
        """
        FFT of array ``fg`` given in G-space.
        """
        ndim, shape = fg.ndim, fg.shape

        if ndim == 1:
            fg = np.reshape(fg, self.shape)
            return self.fft_g2r(fg, fg_ishifted=fg_ishifted).flatten()

        if ndim == 3:
            assert self.size == np.prod(shape[-3:])
            if fg_ishifted: fg = ifftshift(fg)
            fr = ifftn(fg)

        elif ndim > 3:
            assert self.size == np.prod(shape[-3:])
            axes = np.arange(ndim)[-3:]
            if fg_ishifted: fg = ifftshift(fg, axes=axes)
            fr = ifftn(fg, axes=axes)

        else:
            raise NotImplementedError("ndim < 3 are not supported")

        return fr * self.size

def pseudofrosch(x,k,tmax,dt):
    global potential, psi0
    ti=[0]
    fi=[psi0(x)]
    mmax=int((tmax)/dt+1)
    for m in range(1,mmax,1):
        tn=ti[m-1]
        fn=fi[m-1]
        fn=np.array(fn)
        dp1=e**(-1j*potential(x)*dt/2)*fn #
        
        
        ft_fn=ft.fft(dp1) #falschrum
        ft_fn=ft.ifftshift(ft_fn)
        
        dp2=e**(-1j*1/2*k**2*dt)*ft_fn
        
        ift_fn=ft.ifft(dp2)
        ift_fn=ft.ifftshift(ift_fn)
        dp3=e**(-1j*potential(x)*dt/2)*ift_fn
        
        dp3=list(dp3)
        
        ti+=[tn+dt]
        fi=fi+[dp3]
    return [ti,fi]

def gs_mod_gpu(idata,itera=10,osize=256):
    
    
    cut=osize//2
    
    pl=cl.get_platforms()[0]
    devices=pl.get_devices(device_type=cl.device_type.GPU)
    ctx = cl.Context(devices=[devices[0]])
    queue = cl.CommandQueue(ctx)

    plan = Plan(idata.shape, queue=queue,dtype=complex128) #no funciona con "complex128"
    
    src = str(Template(KERNEL).render(
        double_support=all(
            has_double_support(dev) for dev in devices),
        amd_double_support=all(
            has_amd_double_support(dev) for dev in devices)
        ))
    prg = cl.Program(ctx,src).build() 
    

    idata_gpu=cl_array.to_device(queue, ifftshift(idata).astype("complex128"))
    fdata_gpu=cl_array.empty_like(idata_gpu)
    rdata_gpu=cl_array.empty_like(idata_gpu)
    plan.execute(idata_gpu.data,fdata_gpu.data)
    
    mask=exp(2.j*pi*random(idata.shape))
    mask[512-cut:512+cut,512-cut:512+cut]=0
    
    
    idata_gpu=cl_array.to_device(queue, ifftshift(idata+mask).astype("complex128"))
    fdata_gpu=cl_array.empty_like(idata_gpu)
    rdata_gpu=cl_array.empty_like(idata_gpu)
    error_gpu=cl_array.to_device(ctx, queue, zeros(idata_gpu.shape).astype("double"))
    plan.execute(idata_gpu.data,fdata_gpu.data)
    
    e=1000
    ea=1000
    for i in range (itera):
        prg.norm(queue, fdata_gpu.shape, None,fdata_gpu.data)
        plan.execute(fdata_gpu.data,rdata_gpu.data,inverse=True)
        #~ prg.norm1(queue, rdata_gpu.shape,None,rdata_gpu.data,idata_gpu.data,error_gpu.data, int32(cut))
        norm1=prg.norm1
        norm1.set_scalar_arg_dtypes([None, None, None, int32])
        norm1(queue, rdata_gpu.shape,None,rdata_gpu.data,idata_gpu.data,error_gpu.data, int32(cut))
        
        e= sqrt(cl_array.sum(error_gpu).get())/(2*cut)

        #~ if e>ea: 
           #~ 
            #~ break
        #~ ea=e
        plan.execute(rdata_gpu.data,fdata_gpu.data)
    
    fdata=fdata_gpu.get()
    fdata=ifftshift(fdata)
    fdata=exp(1.j*angle(fdata))
    return fdata

def dftups(inp,nor=None,noc=None,usfac=1,roff=0,coff=0):
    """
    Translated from matlab:

     * `Original Source <http://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation/content/html/efficient_subpixel_registration.html>`_
     * Manuel Guizar - Dec 13, 2007
     * Modified from dftus, by J.R. Fienup 7/31/06

    Upsampled DFT by matrix multiplies, can compute an upsampled DFT in just
    a small region.

    This code is intended to provide the same result as if the following
    operations were performed:

      * Embed the array "in" in an array that is usfac times larger in each
        dimension. ifftshift to bring the center of the image to (1,1).
      * Take the FFT of the larger array
      * Extract an [nor, noc] region of the result. Starting with the
        [roff+1 coff+1] element.

    It achieves this result by computing the DFT in the output array without
    the need to zeropad. Much faster and memory efficient than the
    zero-padded FFT approach if [nor noc] are much smaller than [nr*usfac nc*usfac]

    Parameters
    ----------
    usfac : int
        Upsampling factor (default usfac = 1)
    nor,noc : int,int
        Number of pixels in the output upsampled DFT, in units of upsampled
        pixels (default = size(in))
    roff, coff : int, int
        Row and column offsets, allow to shift the output array to a region of
        interest on the DFT (default = 0)
    """
    # this function is translated from matlab, so I'm just going to pretend
    # it is matlab/pylab
    from numpy.fft import ifftshift,fftfreq
    from numpy import pi,newaxis,floor

    nr,nc=np.shape(inp);
    # Set defaults
    if noc is None: noc=nc;
    if nor is None: nor=nr;
    # Compute kernels and obtain DFT by matrix products
    term1c = ( ifftshift(np.arange(nc,dtype='float') - floor(nc/2)).T[:,newaxis] )/nc # fftfreq
    term2c = (( np.arange(noc,dtype='float') - coff  )/usfac)[newaxis,:]              # output points
    kernc=np.exp((-1j*2*pi)*term1c*term2c);

    term1r = ( np.arange(nor,dtype='float').T - roff )[:,newaxis]                # output points
    term2r = ( ifftshift(np.arange(nr,dtype='float')) - floor(nr/2) )[newaxis,:] # fftfreq
    kernr=np.exp((-1j*2*pi/(nr*usfac))*term1r*term2r);
    #kernc=exp((-i*2*pi/(nc*usfac))*( ifftshift([0:nc-1]).' - floor(nc/2) )*( [0:noc-1] - coff ));
    #kernr=exp((-i*2*pi/(nr*usfac))*( [0:nor-1].' - roff )*( ifftshift([0:nr-1]) - floor(nr/2)  ));
    out=np.dot(np.dot(kernr,inp),kernc);
    #return np.roll(np.roll(out,-1,axis=0),-1,axis=1)
    return out 

 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(fft.fftshift(x), y)
     assert_array_almost_equal(fft.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(fft.fftshift(x), y)
     assert_array_almost_equal(fft.ifftshift(y), x)

def plot_sinc_windows():
    fig, ax = plt.subplots(nrows=2, figsize=(8, 8))
    xs = np.linspace(-2, 2., num=200)
    xs4 = np.linspace(-5, 5, num=500)
    y2s = sinc_w(xs, 'lanczos')
    y4s = sinc_w(xs4, 'lanczos', a=5)
    yks = sinc_w(xs4, 'kaiser', a=5, alpha=5)
    yks2 = sinc_w(xs, 'kaiser', a=2, alpha=5)
    ax[0].plot(xs, y2s, "b", label='Lanczos, a=2')
    ax[0].plot(xs4, y4s, 'g', label='Lanczos, a=5')
    ax[0].plot(xs4, yks, "r", label='Kaiser 5, a=5')
    ax[0].plot(xs, yks2, "c", label='Kaiser 5, a=2')
    #ax[0].plot(xs,sinc_w(xs, 'hann'),label='Hann')
    #ax[0].plot(xs,sinc_w(xs, 'kaiser',alpha=5),label='Kaiser=5')
    #ax[0].plot(xs,sinc_w(xs, 'kaiser',alpha=10),label='Kaiser=10')
    #xs4 = np.linspace(-4,4,num=100)
    #ax[0].plot(xs4,sinc_w(xs4, 'lanczos', a = 4), label='Lanczos,a=4')

    ax[0].legend()
    ax[0].set_xlabel(r"$\pm a$")

    #n=400 #zeropadd FFT
    #freqs2 = fftfreq(len(y2s),d=xs[1]-xs[0])
    #freqs4 =fftfreq(400,d=xs4[1]-xs4[0])
    ysh = ifftshift(y2s)
    pady = np.concatenate((ysh[:100], np.zeros((1000,)), ysh[100:]))
    freq2 = fftshift(fftfreq(len(pady), d=0.02))

    ys4h = ifftshift(y4s)

    pad4y = np.concatenate((ys4h[:250], np.zeros((2000,)), ys4h[250:]))

    freq4 = fftshift(fftfreq(len(pad4y), d=0.02))
    fpady = fft(pady)
    fpad4y = fft(pad4y)
    ax[1].plot(freq2, 10 * np.log10(np.abs(fftshift(fpady / fpady[0]))))
    ax[1].plot(freq4, 10 * np.log10(np.abs(fftshift(fpad4y / fpad4y[0]))))

    ysk = ifftshift(yks)
    padk = np.concatenate((ysk[:250], np.zeros((2000,)), ysk[250:]))
    fpadk = fft(padk)
    ax[1].plot(freq4, 10 * np.log10(np.abs(fftshift(fpadk / fpadk[0]))))
    ysk2 = ifftshift(yks2)
    padk2 = np.concatenate((ysk2[:100], np.zeros((1000,)), ysk2[100:]))
    fpadk2 = fft(padk2)
    ax[1].plot(freq2, 10 * np.log10(np.abs(fftshift(fpadk2 / fpadk2[0]))))
    #ax[1].plot(freqs4, fft(ifftshift(
    #ax[1].plot(freqs, get_db_response(xs, 'hann'),label='Hann')
    #ax[1].plot(freqs, get_db_response(xs, 'kaiser',alpha=5),label='Kaiser=5')
    #ax[1].plot(freqs, get_db_response(xs, 'kaiser',alpha=10),label='Kaiser=10')
    #ax[1].legend()
    ax[1].set_ylabel("dB")
    ax[1].set_xlabel("cycles/a")
    plt.show()

def apply_filter(data_array, filter_array,
                 fft_multiplier=None, ifft_multiplier=None,
                 output_multiplier=None, apply_window_func=False,
                 window_shift=None, invert_filter=False):
    """FFT image cube, apply mask and iFFT the output back to image domain.

    Parameters
    ----------
    data_array : ndarray
        Image cube.
    filter_array : ndarray
        Filter to multiply to FFT(data_array). Same shape as data_array.
        Assume (min, max) = (0, 1).
    fft_multiplier : float, optional
        Scalar to multiply to FFT output before masking.
    ifft_multiplier : float, optional
        Scalar to multiply to the masked array before iFFT back.
    output_multiplier : float, optional
        Scalar to multiply to output
    apply_window_func : bool, optional
        Apply blackman-nuttall window function to the first dimension of
        the data_array before FFT.
    window_shift : int
        Shift window function by these many pixels.
        Negative will shift to the left.
    invert_filter: bool, optional
        Invert the filter (1 - filter_array) before applying. Default is False.

    Returns
    -------
    out: ndarray
        Masked image cube.

    """
    if apply_window_func:
        k = nuttall(data_array.shape[0])
        if window_shift:
            k = np.roll(k, window_shift)
        w = (np.ones_like(data_array).T * k).T
    else:
        w = np.ones_like(data_array)
    data_fft = fftshift(fftn(ifftshift(data_array * w)))
    if fft_multiplier is not None:
        data_fft *= fft_multiplier
    if invert_filter:
        filter_array = 1 - filter_array
    data_fft *= filter_array
    out = fftshift(ifftn(ifftshift(data_fft))) / w
    if ifft_multiplier is not None:
        out *= ifft_multiplier
    if output_multiplier is not None:
        out *= output_multiplier
    return out.real

def fast_hessian(img, sigma=3.0):
    """Calculates Hessian in *frequency space* via fft2

    -4π [   μ^2 iL(μ,ν) μν iL(μ,ν)  ]
        [   μν iL(μ,ν)  ν^2 iL(μ,ν) ] 

    where iL(μ,ν) is the inverse 2D FFT of the image

    """
    #ffx = partial(rfft2, s=img.shape)
    #iffx = partial(irfft2, s=img.shape)
    
    ffx = fft2
    iffx = ifft2
    iL = ffx(img)
    
    SHIFT = False

    # coordinate matrices [[ u^2, uv], [uv, v^2]]
    cxx = np.fromfunction(lambda i,j: i**2, iL.shape, dtype='float')
    cxy = np.fromfunction(lambda i,j: i*j, iL.shape, dtype='float')
    cyy = np.fromfunction(lambda i,j: j**2, iL.shape, dtype='float')
    
    if SHIFT:
        cxx = fftshift(cxx)
        cxy = fftshift(cxy)
        cyy = fftshift(cyy)

    # elementwise multiplication
    hxx = -4*np.pi**2 * cxx * iL
    hxy = -4*np.pi**2 * cxy * iL
    hyy = -4*np.pi**2 * cyy * iL
    
    exparg = -(cxx + cyy) / (2*np.pi**2 * sigma**2)
    #A = 1 / (2*np.pi*sigma**2)
    A = 1
    fgauss = A*np.exp(exparg)

    hxx = fgauss * hxx
    hxy = fgauss * hxy
    hyy = fgauss * hyy
    
    Hxx = iffx(hxx)
    Hxy = iffx(hxy)
    Hyy = iffx(hyy)

    if SHIFT:
        Hxx = ifftshift(Hxx)
        Hxy = ifftshift(Hxy)
        Hyy = ifftshift(Hyy)

    return  Hxx, Hxy, Hyy

    def setUp(self):
        # Input
        self.data = hdf5.read_data("example_input.h5")
        C, T, Z, Y, X = self.data[0].shape
        self.data[0].shape = Y, X
        self.small_data = self.data[0][350:400, 325:375]

        # Input drifted by known value
        self.data_drifted = hdf5.read_data("example_drifted.h5")
        C, T, Z, Y, X = self.data_drifted[0].shape
        self.data_drifted[0].shape = Y, X

        # Input drifted by random value
        z = 1j  # imaginary unit
        self.deltar = numpy.random.uniform(-100, 100)
        self.deltac = numpy.random.uniform(-100, 100)
        nr, nc = self.data[0].shape
        array_nr = numpy.arange(-numpy.fix(nr / 2), numpy.ceil(nr / 2))
        array_nc = numpy.arange(-numpy.fix(nc / 2), numpy.ceil(nc / 2))
        Nr = fft.ifftshift(array_nr)
        Nc = fft.ifftshift(array_nc)
        [Nc, Nr] = numpy.meshgrid(Nc, Nr)
        self.data_random_drifted = fft.ifft2(fft.fft2(self.data[0]) * numpy.power(math.e,
        				z * 2 * math.pi * (self.deltar * Nr / nr + self.deltac * Nc / nc)))

        # Noisy inputs
        noise = random.normal(0, 3000, self.data[0].size)
        noise_array = noise.reshape(self.data[0].shape[0], self.data[0].shape[1])

        self.data_noisy = self.data[0] + noise_array
        self.data_drifted_noisy = self.data_drifted[0] + noise_array
        self.data_random_drifted_noisy = self.data_random_drifted + noise_array

        # Small input drifted by random value
        self.small_deltar = numpy.random.uniform(-10, 10)
        self.small_deltac = numpy.random.uniform(-10, 10)
        nr, nc = self.small_data.shape
        array_nr = numpy.arange(-numpy.fix(nr / 2), numpy.ceil(nr / 2))
        array_nc = numpy.arange(-numpy.fix(nc / 2), numpy.ceil(nc / 2))
        Nr = fft.ifftshift(array_nr)
        Nc = fft.ifftshift(array_nc)
        [Nc, Nr] = numpy.meshgrid(Nc, Nr)
        self.small_data_random_drifted = fft.ifft2(fft.fft2(self.small_data) * numpy.power(math.e,
        				z * 2 * math.pi * (self.small_deltar * Nr / nr + self.small_deltac * Nc / nc)))

        # Small noisy inputs
        small_noise = random.normal(0, 3000, self.small_data.size)
        small_noise_array = small_noise.reshape(self.small_data.shape[0], self.small_data.shape[1])

        self.small_data_noisy = self.small_data + small_noise_array
        self.small_data_random_drifted_noisy = self.small_data_random_drifted + small_noise_array

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

def inverse_spectrogram(spec, s_len,
    sample_rate, spec_sample_rate, freq_spacing, min_freq=0, max_freq=None, nstd=6, log=True, noise_level_db=80, rectify=True):
    """turns the complex spectrogram into a signal

    inverts by repeating the process on a string-of-ones
    """

    spec_copy = spec.copy()
    if log:
        spec_copy = 10**(spec_copy)
    spec_tranpose = spec.transpose() # spec_tranpose[time][frequency]

    hnwinlen = len(spec) - 1
    nincrement = int(np.round(float(sample_rate)/spec_sample_rate))

    gauss_t = np.arange(-hnwinlen, hnwinlen+1, 1.0)
    gauss_std = float(2*hnwinlen) / float(nstd)
    gauss_window = np.exp(-gauss_t**2 / (2.0*gauss_std**2)) / (gauss_std*np.sqrt(2*np.pi))
    
    s = np.zeros(s_len + 2*hnwinlen+1)
    w = np.zeros(s_len + 2*hnwinlen+1)

    for i in range(len(spec_tranpose)):
        sample = i * nincrement
        spec_slice = np.concatenate((spec_tranpose[i][:0:-1].conj(), spec_tranpose[i]))
        s[sample:sample+2*hnwinlen+1] += gauss_window * ifft(ifftshift(spec_slice))
        w[sample:sample+2*hnwinlen+1] += gauss_window ** 2
    s /= w
    return s[hnwinlen:hnwinlen+s_len]

def frp_to_fir(frp, fbins=None):
    '''Transform a fringe rate profile to a fir filter.'''
    frp = ifftshift(frp,axes=-1)
    fir = ifft(frp, axis=-1)
    fir = fftshift(fir, axes=-1)
    if fbins is not None: return fir, fftshift(fftfreq(fbins.size, fbins[1] - fbins[0]))
    else: return fir

def delayseq(x, delay_sec:float, fs:int):
    """
    x: input 1-D signal
    delay_sec: amount to shift signal [seconds]
    fs: sampling frequency [Hz]

    xs: time-shifted signal
    """

    assert x.ndim == 1, 'only 1-D signals for now'

    delay_samples = delay_sec*fs
    delay_int = round(delay_samples)

    nfft = nextpow2(x.size+delay_int)

    fbins = 2*pi*ifftshift((arange(nfft)-nfft//2))/nfft

    X = fft(x,nfft)
    Xs = ifft(X*exp(-1j*delay_samples*fbins))

    if isreal(x[0]):
        Xs = Xs.real

    xs = zeros_like(x)
    xs[delay_int:] = Xs[delay_int:x.size]

    return xs

    def compute(self, scene: Scene):
        """ Compute optical irradiance map
        Computation proccedure:
            1) convert radiance to irradiance
            2) apply lens and macular transmittance
            3) apply off-axis fall-off (cos4th)
            4) apply optical transfert function

        Args:
            scene (pyEyeBall.Scene): instance of Scene class, containing the radiance and other scene information

        Examples:
            >>> oi = Optics()
            >>> oi.compute(Scene())
        """
        # set field of view and wavelength samples
        self.fov = scene.fov
        scene.wave = self._wave
        self.dist = scene.dist

        # compute irradiance
        self.photons = pi / (1 + 4 * self.f_number**2 * (1 + abs(self.magnification))**2) * scene.photons

        # apply ocular transmittance
        self.photons *= self.ocular_transmittance

        # apply the relative illuminant (off-axis) fall-off: cos4th function
        x, y = self.spatial_support
        s_factor = np.sqrt(self.image_distance**2 + x**2 + y**2)
        self.photons *= (self.image_distance / s_factor[:, :, None])**4

        # apply optical transfer function of the optics
        for ii in range(self.wave.size):
            otf = fftshift(self.otf(self._wave[ii], self.frequency_support_x, self.frequency_support_y))
            self.photons[:, :, ii] = np.abs(ifftshift(ifft2(otf * fft2(fftshift(self.photons[:, :, ii])))))