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 FFT_Correlation(x,y):
    """
    FFT-based correlation, much faster than numpy autocorr.
    x and y are row-based vectors of arbitrary lengths.
    This is a vectorized implementation of O(N*log(N)) flops.
    """

    lengthx = x.shape[0]
    lengthy = y.shape[0]

    x = np.reshape(x,(1,lengthx))
    y = np.reshape(y,(1,lengthy))

    length = np.array([lengthx, lengthy]).min()
    
    x = x[:length]
    y = y[:length]
    
    fftx = fft(x, 2 * length - 1, axis=1) #pad with zeros
    ffty = fft(y, 2 * length - 1, axis=1)

    corr_xy = fft.ifft(fftx * np.conjugate(ffty), axis=1)
    corr_xy = np.real(fft.fftshift(corr_xy, axes=1)) #should be no imaginary part

    corr_yx = fft.ifft(ffty * np.conjugate(fftx), axis=1)
    corr_yx = np.real(fft.fftshift(corr_yx, axes=1))

    corr = 0.5 * (corr_xy[:,length:] + corr_yx[:,length:]) / range(1,length)[::-1]
    return np.reshape(corr,corr.shape[1])

def B_ell(flat_map, component, bin_size, window=None):
    """
    Return the binned beam function of the given flat map component,
    using ell bins of bin_size. If a window name is given, multiply
    the map component by the window function before taking the 2-D
    DFT.

    The returned beam function is the absolute value of the DFT, so it
    has the same units as the map component. The returned bins array
    contains the left edges of the bins corresponding to the returned
    data: for all i in range(bins.size),
    bins[i] <= data[i] < bins[i+1]
    """
    ell_x, ell_y= np.meshgrid(fftshift(fftfreq(flat_map.x.size, flat_map.dx()/(2*np.pi))),
                              fftshift(fftfreq(flat_map.y.size, flat_map.dy()/(2*np.pi))))
    ell_x = ell_x.T
    ell_y = ell_y.T
    ell_r = np.sqrt(ell_x**2 + ell_y**2)
    ell_bins = np.arange(0, np.max(ell_r), bin_size)
    beam = flat_map[component]
    if window is not None:
        beam *= get_window(window, flat_map.y.size)
        beam *= get_window(window, flat_map.x.size)[:, np.newaxis]
    dft = fftshift(fft2(beam))
    # Shift the indices down by 1 so that they correspond to the data
    # indices. These indices should always be valid indices for the
    # DFT data because r_ell has a lower bound at 0 and max(r_ell)
    # will have index ell_bins.size.
    bin_indices = np.digitize(ell_r.flatten(), ell_bins) - 1
    binned = np.zeros(ell_bins.size) # dtype?
    for i in range(binned.size):
        binned[i] = np.sqrt(np.mean(abs(dft.flatten()[i==bin_indices])**2))
    return ell_bins, binned, dft

def acf2psd(acfall,noiseall,Nr,dns):
    """
    acf all:  Nlag x Nslantrange x real/comp

    """
    assert acfall.ndim in (3,2)

    Nlag = acfall.shape[0]
    spec = empty((Nr,2*Nlag-1),complex128)

    if acfall.ndim == 3: # last dim real,cplx
        acf = (acfall[...,0] + 1j*acfall[...,1]).T / dns / 2.
    elif acfall.ndim == 2 and iscomplex(acfall[0,0]):
        acf = acfall / dns / 2.
    else:
        raise TypeError('is this really ACF? I expect complex 2-D matrix')

    try:
        acf_noise = (noiseall[...,0] + 1j*noiseall[...,1]).T / dns / 2.
    except TypeError:
        acf_noise= None
        spec_noise= 0.
#%% spectrum noise
    if acf_noise is not None:
        spec_noise = zeros(2*Nlag-1,complex128)
        for i in range(Nlag):
            spec_noise += fftshift(fft(append(conj(acf_noise[i,1:][::-1]),acf_noise[i,:])))
        #
        spec_noise = spec_noise / Nlag
#%% spectrum from ACF
    for i in range(Nr):
        spec[i,:] = fftshift(fft(append(conj(acf[i,1:][::-1]), acf[i,:])))-spec_noise

    return spec,acf

def convFFT(x,y):
    nx=len(x)
    xf=np.pad(x,(nx/2,nx/2),mode='constant')
    yf=np.pad(y,(nx/2,nx/2),mode='constant')
    xf=(fft.fft(fft.fftshift(xf)))
    yf=(fft.fft(fft.fftshift(yf)))
    return fft.fftshift(np.real((fft.ifft(xf*yf))))[nx/2:3*nx/2]

    def fft_r2g(self, fr, shift_fg=False):
        """
        FFT of array ``fr`` given in real space.
        """
        ndim, shape = fr.ndim, fr.shape

        if ndim == 1:
            fr = np.reshape(fr, self.shape)
            return self.fft_r2g(fr, shift_fg=shift_fg).flatten()

        elif ndim == 3:
            assert self.size == np.prod(shape[-3:])
            fg = fftn(fr)
            if shift_fg: fg = fftshift(fg)

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

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

        return fg / self.size

 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 _configure(infiles, z, df):
    conf.infiles = infiles
    img_hdr = fits.getheader(conf.infiles[0])
    conf.nx = img_hdr['naxis1']
    conf.ny = img_hdr['naxis2']
    conf.nf = MWA_FREQ_EOR_ALL_80KHZ.size
    conf.dx = np.abs(img_hdr['cdelt1']) * np.pi / 180
    conf.dy = np.abs(img_hdr['cdelt2']) * np.pi / 180
    conf.df = df
    conf.du = 1 / (conf.nx * conf.dx)
    conf.dv = 1 / (conf.ny * conf.dy)
    conf.deta = 1 / (conf.nf * conf.df)
    conf.freq = MWA_FREQ_EOR_ALL_80KHZ
    conf.u = fftshift(fftfreq(conf.nx, conf.dx))
    conf.v = fftshift(fftfreq(conf.ny, conf.dy))
    conf.eta = fftshift(fftfreq(conf.nf, conf.df))
    conf.z = z
    conf.cosmo_d = Cosmo.comoving_transverse_distance(conf.z).value
    conf.cosmo_e = np.sqrt(Cosmo.Om0 * (1 + conf.z) ** 3 + Cosmo.Ok0 *
                           (1 + conf.z) ** 2 + Cosmo.Ode0)
    conf.cosmo_h0 = Cosmo.H0.value * 1e3
    conf.cosmo_c = const.si.c.value
    conf.kx = conf.u * 2 * np.pi / conf.cosmo_d
    conf.ky = conf.v * 2 * np.pi / conf.cosmo_d
    conf.dkx = conf.du * 2 * np.pi / conf.cosmo_d
    conf.dky = conf.dv * 2 * np.pi / conf.cosmo_d
    conf.k_perp = np.sqrt(conf.kx ** 2 + conf.ky[np.newaxis, ...].T ** 2)
    conf.k_par = conf.eta * 2 * np.pi * conf.cosmo_h0 * _F21 * \
        conf.cosmo_e / (conf.cosmo_c * (1 + conf.z) ** 2)

def compacf(acfall,noiseall,Nr,dns,bstride,ti,tInd):
    bstride=bstride.squeeze()
    assert bstride.size==1 #TODO

    Nlag = acfall.shape[2]
    acf  =      zeros((Nr,Nlag),complex64) #NOT empty, note complex64 is single complex float
    spec =      empty((Nr,2*Nlag-1),complex128)
    try:
        acf_noise = zeros((noiseall.shape[3],Nlag),complex64)
        spec_noise= zeros(2*Nlag-1,complex128)
    except AttributeError:
        acf_noise = None
        spec_noise= 0.

    for i in range(tInd[ti]-1,tInd[ti]+1):
        acf += (acfall[i,bstride,:,:,0] + 1j*acfall[i,bstride,:,:,1]).T
        if acf_noise is not None: #must be is not None
            acf_noise += (noiseall[i,bstride,:,:,0] + 1j*noiseall[i,bstride,:,:,1]).T

    acf = acf/dns/(i-(tInd[ti]-1)+1) #NOT /=
    if acf_noise is not None: #must be is not None
        acf_noise = acf_noise/dns / (i-(tInd[ti]-1)+1)
#%% spectrum noise
    if acf_noise is not None:
        for i in range(Nlag):
            spec_noise += fftshift(fft(append(conj(acf_noise[i,1:][::-1]),acf_noise[i,:])))


        spec_noise = spec_noise / Nlag
#%% spectrum from ACF
    for i in range(Nr):
        spec[i,:] = fftshift(fft(append(conj(acf[i,1:][::-1]), acf[i,:])))-spec_noise


    return spec,acf

def sph_fft(x, y):
  dx = x[1] - x[0]
  q = fft.fftshift(fft.fftfreq(len(x), dx)) * 2 * np.pi
  f = -(4*np.pi)*fft.fftshift(np.imag(fft.fft(y*dx * x)))
  f /= q
  f[q == 0] = np.trapz(4*np.pi*y*x**2, x)
  return q, f

def gen_wedge_psf(nx, ny, nf, dx, dy, df, z, out, threads=None):
    u = fftshift(fftfreq(nx, dx * np.pi / 180))
    v = fftshift(fftfreq(ny, dy * np.pi / 180))
    e = fftshift(fftfreq(nf, df))

    E = np.sqrt(Cosmo.Om0 * (1 + z) ** 3 +
                Cosmo.Ok0 * (1 + z) ** 2 + Cosmo.Ode0)
    D = Cosmo.comoving_transverse_distance(z).value
    H0 = Cosmo.H0.value * 1e3
    c = const.c.value
    print(E, D, H0)
    kx = u * 2 * np.pi / D
    ky = v * 2 * np.pi / D
    k_perp = np.sqrt(kx ** 2 + ky[np.newaxis, ...].T ** 2)
    k_par = e * 2 * np.pi * H0 * f21 * E / (c * (1 + z) ** 2)
    arr = np.ones((nf, nx, ny), dtype='complex128')
    for i in range(nf):
        mask = (k_perp > np.abs(k_par[i]) * c * (1 + z) / (H0 * E * D))
        arr[i][mask] = 0
    np.save('kx.npy', kx)
    np.save('ky.npy', ky)
    np.save('kpar.npy', k_par)
    np.save('wedge_window.npy', arr.real)
    fft_arr = fftshift(fftn(ifftshift(arr))).real
    hdu = fits.PrimaryHDU(data=fft_arr)
    hdr_dict = dict(cdelt1=dx, cdelt2=dy, cdelt3=df,
                    crpix1=nx/2, crpix2=ny/2, crpix3=nf/2,
                    crval1=0, crval2=0, crval3=0,
                    ctype1='RA---SIN', ctype2='DEC--SIN', ctype3='FREQ',
                    cunit1='deg', cunit2='deg', cunit3='Hz')
    for k, v in hdr_dict.items():
        hdu.header[k] = v
    hdu.writeto(out, clobber=True)

def create_matching_kernel(source_psf, target_psf, window=None):
    """
    Create a kernel to match 2D point spread functions (PSF) using the
    ratio of Fourier transforms.

    Parameters
    ----------
    source_psf : 2D `~numpy.ndarray`
        The source PSF.  The source PSF should have higher resolution
        (i.e. narrower) than the target PSF.  ``source_psf`` and
        ``target_psf`` must have the same shape and pixel scale.

    target_psf : 2D `~numpy.ndarray`
        The target PSF.  The target PSF should have lower resolution
        (i.e. broader) than the source PSF.  ``source_psf`` and
        ``target_psf`` must have the same shape and pixel scale.

    window : callable, optional
        The window (or taper) function or callable class instance used
        to remove high frequency noise from the PSF matching kernel.
        Some examples include:

        * `~photutils.psf.matching.HanningWindow`
        * `~photutils.psf.matching.TukeyWindow`
        * `~photutils.psf.matching.CosineBellWindow`
        * `~photutils.psf.matching.SplitCosineBellWindow`
        * `~photutils.psf.matching.TopHatWindow`

        For more information on window functions and example usage, see
        :ref:`psf_matching`.

    Returns
    -------
    kernel : 2D `~numpy.ndarray`
        The matching kernel to go from ``source_psf`` to ``target_psf``.
        The output matching kernel is normalized such that it sums to 1.
    """

    # inputs are copied so that they are not changed when normalizing
    source_psf = np.copy(np.asanyarray(source_psf))
    target_psf = np.copy(np.asanyarray(target_psf))

    if source_psf.shape != target_psf.shape:
        raise ValueError('source_psf and target_psf must have the same shape '
                         '(i.e. registered with the same pixel scale).')

    # ensure input PSFs are normalized
    source_psf /= source_psf.sum()
    target_psf /= target_psf.sum()

    source_otf = fftshift(fft2(source_psf))
    target_otf = fftshift(fft2(target_psf))
    ratio = target_otf / source_otf

    # apply a window function in frequency space
    if window is not None:
        ratio *= window(target_psf.shape)

    kernel = np.real(fftshift((ifft2(ifftshift(ratio)))))
    return kernel / kernel.sum()

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 abs_fft(self, points=None, zoom=None,write = 'off'):
		"""
		Fourier transforms the timesignal;
		points is the number of points to transform, if more points given than data points
		the rest is zero padded
	
		absfft(points=4096)
		"""
		realdata = N.array(self.the_result.y[0])
		imdata = N.array(self.the_result.y[1])
		data = realdata + 1j*imdata
		fftdata = F.fftshift(F.fft(data, points))
		absfft = N.sqrt(fftdata.real**2 + fftdata.imag**2)
		# create our x axis
		n = fftdata.size
		self.the_result.x = F.fftshift(F.fftfreq(n, 1.0/self.sampling_rate))
		self.the_result.y[0] = absfft
		self.the_result.y[1] = N.zeros(n)
		if write == 'on':
			return self
		else:
			if zoom is None:
				return self.the_result
			else:
				center, width = zoom
				return self.zoom(self.the_result, center, width)

def delayTransformVisibilities(visList,df,kernel=None,wind='Blackman-Harris'):
    '''
     ARGS:
     visList: nfreqxnvis complex matrix of visibilities
     df: float indicating channel width
     RETURNS:
     ndelayxnvis grid of delay-transformed visibilities
    '''
    nf=visList.shape[0]
    if kernel is None:
        kernel=np.ones(nf)

    nv=visList.shape[1]
    windowGrid=np.zeros_like(visList)
    if wind=='Blackman-Harris':
        window=signal.blackmanharris(nf)
    elif wind=='Nithya':
        window=signal.blackmanharris(nf/2)
        window=signal.fftconvolve(window,window,mode='full')
        window=np.append(window,window[-1])
    window/=window.max()
    window/=np.sqrt(np.mean(np.abs(kernel*window)**2.)) 
    for vNum in range(nv):
        windowGrid[:,vNum]=window*kernel
    return df*fft.fftshift(fft.fft(fft.fftshift(visList*windowGrid,axes=[0]),axis=0),axes=[0])

    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])))))

def fft_2D(R, convert=3e-5, freq_bounds=None):
    """
    Perform a 2D FFT to transform a 3rd order response function defined in the
    rotating frame into a series of 2D spectra.

    First argument must be a MetaArray object wth ticks and rw_freq defined.

    Returns the FFT in another MetaArray object with ticks updated to the
    calculated frequencies (converted using the convert argument which defaults
    to cm to fs).
    """

    # reverses frequency order to keep convention e^{+i \omega t}
    R_2D = fftshift(fft2(ifftshift(R, axes=(0, 2)), axes=(0, 2)),
                    axes=(0, 2))[::-1, :, ::-1]

    dt = [t[1] - t[0] for t in R.ticks]

    freqs = [fftshift(fftfreq(R.shape[axis], dt[axis] * convert))
             + R.rw_freq for axis in (0, 2)]

    if freq_bounds is not None:
        bounds = [bound_indices(ticks, freq_bounds) for ticks in freqs]
        freqs = [freq[bound[0]:bound[1]] for freq, bound in zip(freqs, bounds)]
        R_2D = R_2D[bounds[0][0]:bounds[0][1], :, bounds[1][0]:bounds[1][1]]

    return MetaArray(R_2D, ticks=(freqs[0], R.ticks[1], freqs[1]))

def plot_pixel_effect():
    fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(8, 8))
    #fig = plt.figure()
    #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

    for ax in axs[:, 0]:
        ax.plot(sfreqs, fftshift(tout) / tout[0])
        ax.plot(sfreqs, fftshift(taper))
        ax.plot(sfreqs, 0.0395 * np.sinc(0.0395 * sfreqs))
        ax.plot(sfreqs, 0.0472 * np.sinc(0.0472 * sfreqs))

    for ax in axs[:, 1]:
        ax.plot(sfreqs, 10 * np.log10(np.abs(fftshift(tout) / tout[0])))
        ax.plot(sfreqs, 10 * np.log10(np.abs(fftshift(taper))))
        ax.plot(sfreqs, 10 * np.log10(np.abs(0.0395 * np.sinc(0.0395 * sfreqs))))
        ax.plot(sfreqs, 10 * np.log10(np.abs(0.0472 * np.sinc(0.0472 * sfreqs))))

    axs[0, 0].set_ylabel("Norm amp")
    axs[1, 0].set_ylabel("Norm amp")
    axs[0, 1].set_ylabel("dB")
    axs[1, 1].set_ylabel("dB")
    for ax in axs.flatten():
        ax.set_xlabel(r"cycles/$\lambda$")
    plt.show()

 def plot(self, title, truth=None):
     kwargs = {"interpolation": "nearest",
               "origin": "lower",
               "cmap": "afmhot",
               "vmin": 0.0,
               "vmax": np.max(self.get_real_image())}
     if truth is not None:
         plt.subplot(2,2,3)
         plt.imshow(truth, **kwargs)
         plt.title("truth")
     plt.subplot(2,2,1)
     plt.imshow(self.get_real_image(), **kwargs)
     plt.title("{title}: scores {s1:.1f} {s2:.1f}".format(title=title,
                                                          s1=self.get_score_L1(),
                                                          s2=self.get_score_L2()))
     kwargs["vmin"] = np.log(np.percentile(self.get_data(), 1.))
     kwargs["vmax"] = np.log(np.percentile(self.get_data(), 99.))
     plt.subplot(2,2,4)
     data = np.log(self.get_data().copy())
     data[np.where(self.get_ivar() <= 0)] = kwargs["vmin"]
     plt.imshow(fftshift(data, axes=0), **kwargs)
     plt.title("data")
     plt.subplot(2,2,2)
     plt.title(title)
     plt.imshow(np.log(fftshift(self.get_squared_norm_ft_image(), axes=0)), **kwargs)

def get_spectrum_1d(data_reg,x_reg,y_reg):
    """Compute the 1d power spectrum.
    """
    # remove the mean and squarize
    data_reg-=data_reg.mean()
    jpj,jpi = data_reg.shape
    msize = min(jpj,jpi)
    data_reg = data_reg[:msize-1,:msize-1]
    x_reg = x_reg[:msize-1,:msize-1]
    y_reg = y_reg[:msize-1,:msize-1]
    # wavenumber vector
    x1dreg,y1dreg = x_reg[0,:],y_reg[:,0]
    Ni,Nj = msize-1,msize-1
    dx=npy.int(npy.ceil(x1dreg[1]-x1dreg[0]))
    k_max  = npy.pi / dx
    kx = fft.fftshift(fft.fftfreq(Ni, d=1./(2.*k_max)))
    ky = fft.fftshift(fft.fftfreq(Nj, d=1./(2.*k_max)))
    kkx, kky = npy.meshgrid( ky, kx )
    Kh = npy.sqrt(kkx**2 + kky**2)
    Nmin  = min(Ni,Nj)
    leng  = Nmin/2+1
    kstep = npy.zeros(leng)
    kstep[0] =  k_max / Nmin
    for ind in range(1, leng):
        kstep[ind] = kstep[ind-1] + 2*k_max/Nmin
    norm_factor = 1./( (Nj*Ni)**2 )
    # tukey windowing = tapered cosine window
    cff_tukey = 0.25
    yw=npy.linspace(0, 1, Nj)
    wdw_j = npy.ones(yw.shape)
    xw=npy.linspace(0, 1, Ni)
    wdw_i= npy.ones(xw.shape)
    first_conditioni = xw<cff_tukey/2
    first_conditionj = yw<cff_tukey/2
    wdw_i[first_conditioni] = 0.5 * (1 + npy.cos(2*npy.pi/cff_tukey * (xw[first_conditioni] - cff_tukey/2) ))
    wdw_j[first_conditionj] = 0.5 * (1 + npy.cos(2*npy.pi/cff_tukey * (yw[first_conditionj] - cff_tukey/2) ))
    third_conditioni = xw>=(1 - cff_tukey/2)
    third_conditionj = yw>=(1 - cff_tukey/2)
    wdw_i[third_conditioni] = 0.5 * (1 + npy.cos(2*npy.pi/cff_tukey * (xw[third_conditioni] - 1 + cff_tukey/2)))
    wdw_j[third_conditionj] = 0.5 * (1 + npy.cos(2*npy.pi/cff_tukey * (yw[third_conditionj] - 1 + cff_tukey/2)))
    wdw_ii, wdw_jj = npy.meshgrid(wdw_j, wdw_i, sparse=True)
    wdw = wdw_ii * wdw_jj
    data_reg*=wdw
    #2D spectrum
    cff  = norm_factor
    tempconj=fft.fft2(data_reg).conj()
    tempamp=cff * npy.real(tempconj*fft.fft2(data_reg))
    spec_2d=fft.fftshift(tempamp)
    #1D spectrum
    leng    = len(kstep)
    spec_1d = npy.zeros(leng)
    krange     = Kh <= kstep[0]
    spec_1d[0] = spec_2d[krange].sum()
    for ind in range(1, leng):
        krange = (kstep[ind-1] < Kh) & (Kh <= kstep[ind])
        spec_1d[ind] = spec_2d[krange].sum()
    spec_1d[0] /= kstep[0]
    for ind in range(1, leng):
        spec_1d[ind] /= kstep[ind]-kstep[ind-1]
    return spec_1d, kstep