def deflection_calculation(x, y, topo, rho_t, rho_c, rho_m, Te, E, nu, padding=0):
    """
    Calculates the deflection due to a topographic load for a plate of constant
    thickness Te.
    Uses the equation:
    F[w] = rho_t/(rho_m-rho_c)  phi_e(k)  F[topo]
    """
    ny, nx = np.shape(topo)
    dx = abs(x[0][1] - x[0][0])
    dy = abs(y[1][0] - y[0][0])
    if padding != 0:
        ny_pad, nx_pad = ny*padding, nx*padding
        topo = np.pad(topo, (ny_pad,nx_pad), 'constant', constant_values=0)
    else:
        nx_pad, ny_pad = 0, 0 
        
    fx = fftpack.fftfreq(nx + 2*nx_pad, dx)
    fy = fftpack.fftfreq(ny + 2*ny_pad, dy)
    fx, fy = np.meshgrid(fx, fy)
    k = 2*np.pi*np.sqrt(fx**2 + fy**2)
    
    F_w = rho_t/(rho_m - rho_c)*phi_e(k, Te, rho_c, rho_m, E, nu)*fftpack.fft2(topo)
    w = np.real(fftpack.ifft2(F_w))
    
    if padding != 0:
        w = w[ny_pad:-ny_pad, nx_pad:-nx_pad]
    return w

def make_audio_analysis_plots(infile, prefix='temp', make_plots=True,
                              do_fft=True, fft_sum=None):
    ''' create frequency plot '''
    import numpy as np
    from scipy import fftpack
    from scipy.io import wavfile
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as pl

    if not os.path.exists(infile):
        return -1

    try:
        rate, data = wavfile.read(infile)
    except ValueError:
        print('error reading wav file')
        return -1
    dt_ = 1./rate
    time_ = dt_ * data.shape[0]
    tvec = np.arange(0, time_, dt_)
    sig0 = data[:, 0]
    sig1 = data[:, 1]
    if not tvec.shape == sig0.shape == sig1.shape:
        return -1
    if not do_fft:
        fft_sum_ = float(np.sum(np.abs(sig0)))
        if hasattr(fft_sum, 'value'):
            fft_sum.value = fft_sum_
        return fft_sum_
    if make_plots:
        pl.clf()
        pl.plot(tvec, sig0)
        pl.plot(tvec, sig1)
        xtickarray = range(0, 12, 2)
        pl.xticks(xtickarray, ['%d s' % x for x in xtickarray])
        pl.savefig('%s/%s_time.png' % (HOMEDIR, prefix))
        pl.clf()
    samp_freq0 = fftpack.fftfreq(sig0.size, d=dt_)
    sig_fft0 = fftpack.fft(sig0)
    samp_freq1 = fftpack.fftfreq(sig1.size, d=dt_)
    sig_fft1 = fftpack.fft(sig1)
    if make_plots:
        pl.clf()
        pl.plot(np.log(np.abs(samp_freq0)+1e-9), np.abs(sig_fft0))
        pl.plot(np.log(np.abs(samp_freq1)+1e-9), np.abs(sig_fft1))
        pl.xlim(np.log(10), np.log(40e3))
        xtickarray = np.log(np.array([20, 1e2, 3e2, 1e3, 3e3, 10e3, 30e3]))
        pl.xticks(xtickarray, ['20Hz', '100Hz', '300Hz', '1kHz', '3kHz',
                               '10kHz', '30kHz'])
        pl.savefig('%s/%s_fft.png' % (HOMEDIR, prefix))
        pl.clf()

        run_command('mv %s/%s_time.png %s/%s_fft.png %s/public_html/videos/'
                    % (HOMEDIR, prefix, HOMEDIR, prefix, HOMEDIR))

    fft_sum_ = float(np.sum(np.abs(sig_fft0)))
    if hasattr(fft_sum, 'value'):
        fft_sum.value = fft_sum_
    return fft_sum_

 def test_definition(self):
     x = [0,1,2,3,4,-4,-3,-2,-1]
     assert_array_almost_equal(9*fftfreq(9),x)
     assert_array_almost_equal(9*pi*fftfreq(9,pi),x)
     x = [0,1,2,3,4,-5,-4,-3,-2,-1]
     assert_array_almost_equal(10*fftfreq(10),x)
     assert_array_almost_equal(10*pi*fftfreq(10,pi),x)

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 invert_vort(uc, dx=dx, dy=dy, nx=nx, ny=ny, geom=tad.geom):

    ucv = geom.validview(uc)
    vort = ucv[0]
    u = ucv[1]
    v = ucv[2]

    f = fft2(vort)

    nx, ny = vort.shape

    scal_y = 2*pi/dy/ny
    scal_x = 2*pi/dx/nx

    k = fftfreq(nx, 1/nx)[:,None] * 1j * scal_x
    l = fftfreq(ny, 1/ny)[None,:] * 1j * scal_y


    lapl = k**2 + l**2
    lapl[0,0] = 1.0

    psi = f/lapl
    u[:] = -real(ifft2(psi * l))
    v[:] = real(ifft2(psi * k))

    return uc

    def icwt2d(self, da=0.25):
        '''
        Inverse bi-dimensional continuous wavelet transform as in Wang and
        Lu (2010), equation [5].

        Parameters
        ----------
        da : float, optional
            Spacing in the frequency axis.
        '''
        if self.Wf is None:
            raise TypeError("Run cwt2D before icwt2D")
        m0, l0, k0 = self.Wf.shape

        if m0 != self.scales.size:
            raise Warning('Scale parameter array shape does not match\
                           wavelet transform array shape.')
        # Calculates the zonal and meridional wave numters.
        L, K = 2 ** int(np.ceil(np.log2(l0))), 2 ** int(np.ceil(np.log2(k0)))
        # Calculates the zonal and meridional wave numbers.
        l, k = fftfreq(L, self.dy), fftfreq(K, self.dx)
        # Creates empty inverse wavelet transform array and fills it for every
        # discrete scale using the convolution theorem.
        self.iWf = np.zeros((m0, L, K), 'complex')
        for i, an in enumerate(self.scales):
            psi_ft_bar = an * self.wavelet.psi_ft(an * k, an * l)
            W_ft = fftn(self.Wf[i, :, :], s=(L, K))
            self.iWf[i, :, :] = ifftn(W_ft * psi_ft_bar, s=(L, K)) *\
                da / an ** 2.

        self.iWf = self.iWf[:, :l0, :k0].real.sum(axis=0) / self.wavelet.cpsi

        return self

    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 __init__(self, cube, header, phys_units=False):
        super(VCS, self).__init__()

        self.header = header
        self.cube = cube
        self.fftcube = None
        self.correlated_cube = None
        self.ps1D = None
        self.phys_units = phys_units

        if np.isnan(self.cube).any():
            self.cube[np.isnan(self.cube)] = 0
            # Feel like this should be more specific
            self.good_channel_count = np.sum(self.cube[0, :, :] != 0)
        else:
            self.good_channel_count = float(
                self.cube.shape[1] * self.cube.shape[2])

        # Lazy check to make sure we have units of km/s
        if np.abs(self.header["CDELT3"]) > 1:
            self.vel_to_pix = np.abs(self.header["CDELT3"]) / 1000.
        else:
            self.vel_to_pix = np.abs(self.header["CDELT3"])

        self.vel_channels = np.arange(1, self.cube.shape[0], 1)

        if self.phys_units:
            self.vel_freqs = np.abs(
                fftfreq(self.cube.shape[0])) / self.vel_to_pix
        else:
            self.vel_freqs = np.abs(fftfreq(self.cube.shape[0]))

def direct_shift(x,a,period=None):
    n = len(x)
    if period is None:
        k = fftfreq(n)*1j*n
    else:
        k = fftfreq(n)*2j*pi/period*n
    return ifft(fft(x)*exp(k*a)).real

    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 analyze_audio(prefix):
    output = []
    rate, data = wavfile.read('%s.wav' % prefix)
    dt = 1./rate
    T = dt * data.shape[0]
    output.append('%s %s' % (dt, T))

    tvec = np.arange(0, T, dt)
    sig0 = data[:, 0]
    sig1 = data[:, 1]

    output.append('%s %s' % (np.sum(sig0), np.sum(sig1)))

    plt.clf()
    plt.plot(tvec, sig0)
    plt.plot(tvec, sig1)
    xtickarray = range(0, 12, 2)
    plt.xticks(xtickarray, ['%d s' % x for x in xtickarray])
    plt.savefig('%s_time.png' % prefix)

    plt.clf()
    samp_freq0 = fftpack.fftfreq(sig0.size, d=dt)
    sig_fft0 = fftpack.fft(sig0)
    samp_freq1 = fftpack.fftfreq(sig1.size, d=dt)
    sig_fft1 = fftpack.fft(sig1)
    plt.plot(np.log(np.abs(samp_freq0)+1e-9), np.abs(sig_fft0))
    plt.plot(np.log(np.abs(samp_freq1)+1e-9), np.abs(sig_fft1))
    plt.xlim(np.log(10), np.log(40e3))
    xtickarray = np.log(np.array([20, 1e2, 3e2, 1e3, 3e3, 10e3, 30e3]))
    plt.xticks(xtickarray, ['20Hz', '100Hz', '300Hz', '1kHz', '3kHz', '10kHz',
                            '30kHz'])
    plt.savefig('%s_freq.png' % prefix)

    return '\n'.join(output)

def get_mps(t, freq, spec):
    "Computes the MPS of a spectrogram (idealy a log-spectrogram) or other REAL time-freq representation"
    mps = fftshift(fft2(spec))
    amps = np.real(mps * np.conj(mps))
    nf = mps.shape[0]
    nt = mps.shape[1]
    wfreq = fftshift(fftfreq(nf, d=freq[1] - freq[0]))
    wt = fftshift(fftfreq(nt, d=t[1] - t[0]))
    return wt, wfreq, mps, amps

 def mkgrad(self,data):
   """Calculate gradient by convolution with appropriate gaussian filter"""
   n1=data.shape[0]
   if(not self.n or n1 != self.n):
     self.n= n1
     self.grad_multiplyer= fftfreq(n1)*numpy.exp(-4*fftfreq(n1)**2)
   data_f = fft(data)
   grad_f = data_f*self.grad_multiplyer
   return numpy.absolute(ifft(grad))

 def makeImageFromSf(self, sfx, sf, reduceNoise=True):
     """Invert SF all the way back to ImageI."""
     self.invertSf(sfx, sf)
     self.invertAcovf1d(reduceNoise=reduceNoise)
     self.invertAcovf2d(useI=True)
     self.invertPsd2d(useI=True)
     self.invertFft(useI=True)
     self.xfreq = fftpack.fftfreq(self.nx, 1.0)
     self.yfreq = fftpack.fftfreq(self.ny, 1.0)
     return

    def test_Eady(self, atol=5e-2, nz=20, Ah=0.):
        """ Eady setup
        """
        ###########
        # prepare parameters for Eady
        ###########
        nz = nz
        zin = np.arange(nz+1, dtype=np.float64)/nz
        N2 = np.full(nz, 1.)
        f0 = 1.
        beta = 0.
        Nx = 10
        Ny = 1
        dx = 1e-1
        dy = 1e-1
        k = fft.fftshift( fft.fftfreq(Nx, dx) )
        l = fft.fftshift( fft.fftfreq(Ny, dy) )
        vbar = np.zeros(nz+1)
        ubar = zin
        etax = np.zeros(2)
        etay = np.zeros(2)

        z, growth_rate, vertical_modes = \
            modes.instability_analysis_from_N2_profile(
                .5*(zin[1:]+zin[:-1]), N2, f0, beta, k, l,
                zin, ubar, vbar, etax, etay, depth=1., sort='LI', num=2
        )

        self.assertEqual(nz+1, vertical_modes.shape[0],
            msg='modes array must be in the right shape')

        self.assertTrue(np.all( np.diff(
                    growth_rate.reshape((growth_rate.shape[0], len(k)*len(l))).imag.max(axis=1) ) <= 0.),
            msg='imaginary part of modes should be descending')

        mode_amplitude1 = (np.absolute(vertical_modes[:, 0])**2).sum(axis=0)
        self.assertTrue(np.allclose(1., mode_amplitude1),
            msg='mode1 should be normalized to amplitude of 1 at all horizontal wavenumber points')

        mode_amplitude2 = (np.absolute(vertical_modes[:, 1])**2).sum(axis=0)
        self.assertTrue(np.allclose(1., mode_amplitude2),
            msg='mode2 should be normalized to amplitude of 1 at all horizontal wavenumber points')

        #########
        # Analytical solution for Eady growth rate
        #########
        growthEady = np.zeros(len(k))
        for i in range(len(k)):
            if (k[i]==0) or ((np.tanh(.5*k[i])**-1 - .5*k[i]) * (.5*k[i] - np.tanh(.5*k[i])) < 0):
                growthEady[i] = 0.
            else:
                growthEady[i] = ubar.max() * np.sqrt( (np.tanh(.5*k[i])**-1 - .5*k[i]) * (.5*k[i] - np.tanh(.5*k[i])) )

        self.assertTrue( np.allclose(growth_rate.imag[0, 0, :], growthEady, atol=atol),
            msg='The numerical growth rates should be close to the analytical Eady solution' )

    def __init__(self, datain, spd, tim_taper=0.1):
        """Arguments:
        
       'datain'    -- the data to be filtered. dimension must be (time, lat, lon)

       'spd'       -- samples per day

       'tim_taper' -- tapering ratio by cos. applay tapering first and last tim_taper%
                      samples. default is cos20 tapering

                      """
        ntim, nlat, nlon = datain.shape

        #remove the lowest three harmonics of the seasonal cycle (WK99, WKW03)
##         if ntim > 365*spd/3:
##             rf = fftpack.rfft(datain,axis=0)
##             freq = fftpack.rfftfreq(ntim*spd, d=1./float(spd))
##             rf[(freq <= 3./365) & (freq >=1./365),:,:] = 0.0     #freq<=3./365 only??
##             datain = fftpack.irfft(rf,axis=0)

        #remove dominal trend
        data = signal.detrend(datain, axis=0)

        #tapering
        if tim_taper == 'hann':
            window = signal.hann(ntim)
            data = data * window[:,NA,NA]
        elif tim_taper > 0:
        #taper by cos tapering same dtype as input array
            tp = int(ntim*tim_taper)
            window = numpy.ones(ntim, dtype=datain.dtype)
            x = numpy.arange(tp)
            window[:tp] = 0.5*(1.0-numpy.cos(x*pi/tp))
            window[-tp:] = 0.5*(1.0-numpy.cos(x[::-1]*pi/tp))
            data = data * window[:,NA,NA]

        #FFT
        self.fftdata = fftpack.fft2(data, axes=(0,2))

        #Note
        # fft is defined by exp(-ikx), so to adjust exp(ikx) multipried minus         
        wavenumber = -fftpack.fftfreq(nlon)*nlon
        frequency = fftpack.fftfreq(ntim, d=1./float(spd))
        knum, freq = numpy.meshgrid(wavenumber, frequency)

        #make f<0 domain same as f>0 domain
        #CAUTION: wave definition is exp(i(k*x-omega*t)) but FFT definition exp(-ikx)
        #so cahnge sign
        knum[freq<0] = -knum[freq<0]
        freq = numpy.abs(freq)
        self.knum = knum
        self.freq = freq

        self.wavenumber = wavenumber
        self.frequency = frequency

def mps(spectrogram, df, dt):
    """
        Compute the modulation power spectrum for a given spectrogram.
    """

    #normalize and mean center the spectrogram
    sdata = copy.copy(spectrogram)
    sdata /= sdata.max()
    sdata -= sdata.mean()

    #take the 2D FFT and center it
    smps = fft2(sdata)
    smps = fftshift(smps)

    #compute the log amplitude
    mps_logamp = 20*np.log10(np.abs(smps)**2)
    mps_logamp[mps_logamp < 0.0] = 0.0

    #compute the phase
    mps_phase = np.angle(smps)

    #compute the axes
    nf = mps_logamp.shape[0]
    nt = mps_logamp.shape[1]
    spectral_freq = fftshift(fftfreq(nf, d=df))
    temporal_freq = fftshift(fftfreq(nt, d=dt))

    """
    nb = sdata.shape[1]
    dwf = np.zeros(nb)
    for ib in range(int(np.ceil((nb+1)/2.0))+1):
        posindx = ib
        negindx = nb-ib+2
        print 'ib=%d, posindx=%d, negindx=%d' % (ib, posindx, negindx)
        dwf[ib]= (ib-1)*(1.0/(df*nb))
        if ib > 1:
            dwf[negindx] =- dwf[ib]

    nt = sdata.shape[0]
    dwt = np.zeros(nt)
    for it in range(0, int(np.ceil((nt+1)/2.0))+1):
        posindx = it
        negindx = nt-it+2
        print 'it=%d, posindx=%d, negindx=%d' % (it, posindx, negindx)
        dwt[it] = (it-1)*(1.0/(nt*dt))
        if it > 1 :
            dwt[negindx] = -dwt[it]

    spectral_freq = dwf
    temporal_freq = dwt
    """

    return temporal_freq,spectral_freq,mps_logamp,mps_phase

    def __init__(self,shape,lengths):
        self.lengths=lengths
        self.shape=shape

        self.Inv_Laplacian=(np.dstack(np.meshgrid(fftpack.fftfreq(self.shape[1],1.0/self.shape[1]),
                                       fftpack.fftfreq(self.shape[0],1.0/self.shape[0])))**2).sum(-1)
        self.Inv_Laplacian[self.Inv_Laplacian<1.0]=1.0
        self.Inv_Laplacian=1.0
        self.Inv_Laplacian/=self.lengths.area_lat_lon

        #Define the vector calculus in spherical harmonics:
        self.spharm=horizontal_vector_calculus_spherical_harmonics(self.shape,self.lengths)
        return

    def solve_sqg(self):
        import scipy.fftpack as fft
        import numpy as np
        from math import pi
        
        self.precondition()
        
        dx = self.dx
        dy = self.dy
        rhos = self.ssd
        
        bhat = fft.fft2( - 9.81 * rhos / self.rho0)  # calculate bouyance 
        ny, nx = rhos.shape
        nz = self.nz
        k = 2 * pi * fft.fftfreq(nx)
        l = 2 * pi * fft.fftfreq(ny)

        ipsihat = np.zeros((nz+3, ny, nx))*complex(0, 0)
        
        Q = np.zeros((nz + 1, 1), dtype='float64'); Q[[0,-1]] = 0.0 # for interior PV, no used in this version
        
        # cutoff value
        ck, cl = 2 * pi / self.filterL, 2 * pi / self.filterL
        # loop through wavenumbers
        bhats = np.zeros_like(bhat)
        for ik in np.arange(k.size):
            for il in np.arange(l.size):
                wv2 = ((k[ik] / dx[il, 0]) ** 2 + (l[il] / dy[0, ik]) ** 2)
                if wv2 > (ck * ck + cl * cl):
                    bhats[il,ik] = bhat[il,ik]
                    right = - bhat[il, ik] / self.f0 * self.Rp
                    left = self.M - wv2 * np.eye(self.nz+1)
                    ipsihat[1:-1, il, ik] = np.linalg.solve(left, right).flatten()
                else:
                    print 'skip k(ik,il)', ik, il, "wv2 = ", wv2
        
        for k in range(1,nz+2):
            ipsihat[k, :, :] = (fft.ifft2(ipsihat[k, :, :]))
        
        if self.bottomboundary == 'psi=0':
            self.psis = np.r_[(np.real(ipsihat)), np.zeros((1,ny,nx))]
        else:
            self.psis = np.real(ipsihat)
            self.psis[0,:,:]= self.psis[1,:,:]
            self.psis[-1,:,:]=self.psis[-2,:,:]-self.dzc[-1]*np.real(fft.ifft2(-bhats))/self.f0
            
        self.rhos = self.psi2rho(self.psis)
        self.us, self.vs = psi2uv(self.lon, self.lat, self.psis)
        
        return

    def __init__(self, shape, space):
        "docstring"
        nx, ny = shape
        dx, dy = space

        scal_y = 2 * pi / dy / ny
        scal_x = 2 * pi / dx / nx

        k = fftfreq(nx, 1 / nx)[:, None] * 1j * scal_x
        l = fftfreq(ny, 1 / ny)[None, :] * 1j * scal_y

        lapl = k**2 + l**2
        lapl[0, 0] = 1.0

        self.k, self.l, self.lapl = k, l, lapl