def flux_qg(q, parms):

    # - (u + U) q_x - (q_y + Q_y) v
    qe = np.vstack((q,-np.flipud(q)))
    qe_hat = fftn(qe)

    # Compute gradient of PV
    q_x = (ifftn( parms.ikx*qe_hat)).real
    q_y = (ifftn( parms.iky*qe_hat)).real

    # Compute streamfunction
    psie_hat = parms.K2Fi*qe_hat
    psi = (ifftn(psie_hat)).real

    # Compute physical velocities
    u = (ifftn(-parms.iky*psie_hat)).real
    v = (ifftn( parms.ikx*psie_hat)).real

    # Restrict to physical domain
    q_x = q_x[0:parms.Ny,:]
    q_y = q_y[0:parms.Ny,:]
    u   = u[0:parms.Ny,:]
    v   = v[0:parms.Ny,:]
    psi = psi[0:parms.Ny,:]

    # Compute flux
    flux = - (u + parms.U)*q_x - (q_y + parms.Q_y)*v

    # FJP: energy should include potential energy
    energy = 0.5*np.mean(u**2 + v**2) + np.mean(parms.F*psi**2)
    enstr  = np.mean(q**2)
    mass   = np.mean(psi)

    return flux, energy, enstr, mass

 def daughter(self, scale, N, dt=1):
     '''
     Returns a daughter wavelet to be multiplied with the
     Fourier-transformed time series.
     
     Parameters
     ----------
     scale : float
         Scale of the wavelet.
     N : int
         Number of observations in the series being transformed.
     dt : int
         Number of observations per unit time.
         
     Returns
     -------
     daughter : ndarray
         "Daughter" wavelet (the Fourier transform of the
         mother wavelet of the appropriate scale and length)
     '''
     k = sp.arange(int(N/2)) * 2 * sp.pi / (N * dt)
     k = sp.hstack((k, -sp.flipud(k)))
     if len(k) == N + 1:
         k = k[1: ]
     expnt = -(scale * k - self.k0)**2 / 2. * (k > 0)
     norm = sp.sqrt(scale * k[1]) * sp.pi**(-0.25) * sp.sqrt(N) # total energy = N (Eqn. 7)
     daughter = norm * sp.exp(expnt)
     daughter = daughter * (k >0)
     return daughter   

def prob4(filename='saw.wav', new_rate = 11025, outfile='prob4.wav'):
    """Down-samples a given .wav file to a new rate and saves the resulting
    signal as another .wav file.
    
    Parameters
    ----------
    filename : string, optional
        The name of the .wav sound file to be down-sampled.
        Defaults to 'saw.wav'.
    new_rate : integer, optional
        The down-sampled rate. Defaults to 11025.
    outfile : string, optional
        The name of the new file. Defaults to prob4.wav.

    Returns
    -------
    None
    """
    old_rate, in_sig = wavfile.read(filename)
    fin = fftw.fft(sp.float32(in_sig))
    # Use if scipy_fftpack is unavailable
    # fin = sp.fft(sp.float32(in_sig))
    nsiz = sp.floor(in_sig.size * new_rate / old_rate)
    nsizh = sp.floor(nsiz / 2)
    fout = sp.zeros(nsiz) + 0j
    fout[0:nsizh] = fin[0:nsizh]
    fout[nsiz-nsizh+1:] = sp.conj(sp.flipud(fout[1:nsizh]))
    out = sp.real(sp.ifft(fout))
    out = sp.int16(out/sp.absolute(out).max() * 32767)
    plot_signal(filename)
    wavfile.write('prob4.wav',new_rate,out)
    print ""; plot_signal('prob4.wav')

def down_sample(filename, new_rate, outputfile=None):
    """
    Create a down-sampled copy of the provided .wav file.  Unless overridden, the output
        file will be of the form "down_<orginalname>.wav"
        
    Parameters
    ----------
    filename : string
        input .wav file
    new_rate : int
        sample rate of output file
    outputfile : string
        name of output file
    """

    if outputfile is None:
        outputfile = "down_" + filename

    old_rate, in_sig = wavfile.read(filename)
    in_sig = sp.float32(in_sig)
    fin = sp.fft(in_sig)
    nsiz = sp.floor(in_sig.size * new_rate / old_rate)
    nsizh = sp.floor(nsiz / 2)
    fout = sp.zeros(nsiz)
    fout = fout + 0j
    fout[0:nsizh] = fin[0:nsizh]
    fout[nsiz - nsizh + 1 :] = sp.conj(sp.flipud(fout[1:nsizh]))
    out = sp.ifft(fout)
    out = sp.real(out)  # Take the real component of the signal
    out = sp.int16(out / sp.absolute(out).max() * 32767)
    wavfile.write(outputfile, new_rate, out)

def getFluxes(val_mat, direction_mat, dist_mat, duxdy_mat, out_flux, inc):

	import scipy;
	import math;

	speed_factor   = 1;
	angle_factor   = 1;
	inc_factor     = 1;
	dist_factor    = 1;
	strain_factor  = 1;

	duxdy_mat = duxdy_mat / (sum(duxdy_mat[~scipy.isnan(duxdy_mat)]));

	cell_angles  = scipy.flipud(scipy.array([[-1 * math.pi / 4, -1 * math.pi / 2, -3 * math.pi / 4], [0, scipy.nan, math.pi], [math.pi / 4, math.pi / 2, 3 * math.pi / 4]]));
#	cell_angles  = scipy.flipud(scipy.array([[3 * math.pi / 4, 1 * math.pi / 2, 1 * math.pi / 4], [math.pi, scipy.nan, 0], [-3 * math.pi / 4, -1 * math.pi / 2, -1 * math.pi / 4]]));

	cell_incs    = scipy.array([[(inc**2 + inc**2)**0.5, inc, (inc**2 + inc**2)**0.5], [inc, scipy.nan, inc], [(inc**2 + inc**2)**0.5, inc, (inc**2 + inc**2)**0.5]]);
	cell_incs    = (1 / cell_incs**inc_factor);
	cell_incs    = cell_incs / sum(cell_incs[~scipy.isnan(cell_incs)]);

	vels_in      = scipy.cos(cell_angles - direction_mat);
	vels_in[1,1] = scipy.nan;
	vels_in[vels_in < 0.00001] = scipy.nan;
	vels_in      = vels_in**angle_factor * val_mat**speed_factor * dist_mat**dist_factor * (1 / duxdy_mat**strain_factor) * cell_incs;
	in_fluxes    = (vels_in / sum(vels_in[~scipy.isnan(vels_in)]) * out_flux);

        return in_fluxes;

 def readData(self,ob):
     f = open(ob.filename,'rb')
     f.seek(int(self._header['STM image list']['Data offset']))
     data = f.read(int(self._header['STM image list']['Data length']))
     ob.d = scipy.fromstring(data,dtype=scipy.int16)
     ob.d.shape = ob.XRes, ob.YRes
     ob.d = scipy.flipud(ob.d)

def KramersKronigFFT(ImX_A):
	'''	Hilbert transform used to calculate real part of a function from its imaginary part
	uses piecewise cubic interpolated integral kernel of the Hilbert transform
	use only if len(ImX_A)=2**m-1, uses fft from scipy.fftpack  '''
	X_A = sp.copy(ImX_A)
	N = int(len(X_A))
	## be careful with the data type, orherwise it fails for large N
	if N > 3e6: A = sp.arange(3,N+1,dtype='float64')
	else:       A = sp.arange(3,N+1)  
	X1 = 4.0*sp.log(1.5)
	X2 = 10.0*sp.log(4.0/3.0)-6.0*sp.log(1.5)
	## filling the kernel
	if N > 3e6: Kernel_A = sp.zeros(N-2,dtype='float64')
	else:       Kernel_A = sp.zeros(N-2)
	Kernel_A = (1-A**2)*((A-2)*sp.arctanh(1.0/(1-2*A))+(A+2)*sp.arctanh(1.0/(1+2*A)))\
	+((A**3-6*A**2+11*A-6)*sp.arctanh(1.0/(3-2*A))+(A+3)*(A**2+3*A+2)*sp.arctanh(1.0/(2*A+3)))/3.0
	Kernel_A = sp.concatenate([-sp.flipud(Kernel_A),sp.array([-X2,-X1,0.0,X1,X2]),Kernel_A])/sp.pi
	## zero-padding the functions for fft
	ImXExt_A = sp.concatenate([X_A[int((N-1)/2):],sp.zeros(N+2),X_A[:int((N-1)/2)]])
	KernelExt_A = sp.concatenate([Kernel_A[N:],sp.zeros(1),Kernel_A[:N]])
	## performing the fft
	ftReXExt_A = -fft(ImXExt_A)*fft(KernelExt_A)
	ReXExt_A = sp.real(ifft(ftReXExt_A))
	ReX_A = sp.concatenate([ReXExt_A[int((3*N+3)/2+1):],ReXExt_A[:int((N-1)/2+1)]])
	return ReX_A

def XIntegralsFFT(GF_A,Bubble_A,Lambda,BubZero):
	''' calculate X integral to susceptibilities using FFT '''
	N = int((len(En_A)-1)/2)
	Kappa_A  = TwoParticleBubble(GF_A,GF_A**2,'eh')
	Bubble_A = TwoParticleBubble(GF_A,GF_A,'eh')
	#print(Kappa_A[N],Bubble_A[N])
	V_A   = 1.0/(1.0+Lambda*Bubble_A)
	KV_A  = Lambda*Kappa_A*V_A**2
	KmV_A = Lambda*sp.flipud(sp.conj(Kappa_A))*V_A**2
	## zero-padding the arrays
	exFD_A  = sp.concatenate([FD_A[N:],sp.zeros(2*N+2),FD_A[:N+1]])
	ImGF_A  = sp.concatenate([sp.imag(GF_A[N:]),sp.zeros(2*N+2),sp.imag(GF_A[:N+1])])
	ImGF2_A = sp.concatenate([sp.imag(GF_A[N:]**2),sp.zeros(2*N+2),sp.imag(GF_A[:N+1]**2)])
	ImV_A   = sp.concatenate([sp.imag(V_A[N:]),sp.zeros(2*N+2),sp.imag(V_A[:N+1])])
	ImKV_A  = sp.concatenate([sp.imag(KV_A[N:]),sp.zeros(2*N+2),sp.imag(KV_A[:N+1])])
	ImKmV_A = sp.concatenate([sp.imag(KmV_A[N:]),sp.zeros(2*N+2),sp.imag(KmV_A[:N+1])])
	## performing the convolution
	ftImX11_A = -sp.conj(fft(exFD_A*ImV_A))*fft(ImGF2_A)*dE
	ftImX12_A =  fft(exFD_A*ImGF2_A)*sp.conj(fft(ImV_A))*dE
	ftImX21_A = -sp.conj(fft(exFD_A*ImKV_A))*fft(ImGF_A)*dE
	ftImX22_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKV_A))*dE
	ftImX31_A = -sp.conj(fft(exFD_A*ImKmV_A))*fft(ImGF_A)*dE
	ftImX32_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKmV_A))*dE
	## inverse transform
	ImX1_A =  sp.real(ifft(ftImX11_A+ftImX12_A))/sp.pi
	ImX2_A =  sp.real(ifft(ftImX21_A+ftImX22_A))/sp.pi
	ImX3_A = -sp.real(ifft(ftImX31_A+ftImX32_A))/sp.pi
	ImX1_A =  sp.concatenate([ImX1_A[3*N+4:],ImX1_A[:N+1]])
	ImX2_A =  sp.concatenate([ImX2_A[3*N+4:],ImX2_A[:N+1]])
	ImX3_A =  sp.concatenate([ImX3_A[3*N+4:],ImX3_A[:N+1]])
	## getting real part from imaginary
	X1_A = KramersKronigFFT(ImX1_A) + 1.0j*ImX1_A + BubZero # constant part !!!
	X2_A = KramersKronigFFT(ImX2_A) + 1.0j*ImX2_A
	X3_A = KramersKronigFFT(ImX3_A) + 1.0j*ImX3_A
	return [X1_A,X2_A,X3_A]

 def bulk_bands_calculator(self,s,sub,kx,ky,kz):
   ''' Calculate the band energies for the specified kx, ky, and kz values.
       The 3x3 Hamiltonian for wurtzite crystals is used for the valence,
       while a 1x1 Hamiltonian is used for the conduction band. The model is
       from the chapter by Vurgaftman and Meyer in the book by Piprek. 
   '''
   E = scipy.zeros((4,len(s.Eg0)))   
   E[0,:] = s.Eg0+s.delcr+s.delso/3+\
               hbar**2/(2*s.mepara)*(kx**2+ky**2)+\
               hbar**2/(2*s.meperp)*(kz**2)+\
               (s.a1+s.D1)*s.epszz+(s.a2+s.D2)*(s.epsxx+s.epsyy)
   L = hbar**2/(2*m0)*(s.A1*kz**2+s.A2*(kx+ky)**2)+\
       s.D1*s.epszz+s.D2*(s.epsxx+s.epsyy)
   T = hbar**2/(2*m0)*(s.A3*kz**2+s.A4*(kx+ky)**2)+\
       s.D3*s.epszz+s.D4*(s.epsxx+s.epsyy)
   F = s.delcr+s.delso/3+L+T
   G = s.delcr-s.delso/3+L+T
   K = hbar**2/(2*m0)*s.A5*(kx+1j*ky)**2+s.D5*(s.epsxx-s.epsyy)
   H = hbar**2/(2*m0)*s.A6*(kx+1j*ky)*kz+s.D6*(s.epsxz)
   d = scipy.sqrt(2)*s.delso/3
   for ii in range(len(s.Eg0)):
     mat = scipy.matrix([[    F[ii],     K[ii],       -1j*H[ii]      ],
                         [    K[ii],     G[ii],       -1j*H[ii]+d[ii]],
                         [-1j*H[ii], -1j*H[ii]+d[ii],     L[ii]      ]])
     w,v = scipy.linalg.eig(mat)
     E[1:,ii] = scipy.flipud(scipy.sort(scipy.real(w)))
   return E

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 plot_patch(ls,H,sigma68,sigma95,sigma99):
    print("Plotting confidence intervals")
    x = sp.append(ls, sp.flipud(ls))
    
    # Initiate figure
#    plt.clf()
#    fig = plt.figure(1)
    
    # Columnwise mean of H:
    Hmean = H.mean(axis=0)

    y68_l = scale(ls,-sigma68[:,0]+Hmean)
    y68_u = scale(ls,sigma68[:,1]+Hmean)
    y68 = sp.append(y68_l,sp.flipud(y68_u))
    
    y95_l = scale(ls,-sigma95[:,0]+Hmean)
    y95_u = scale(ls,sigma95[:,1]+Hmean)
    y95 = sp.append(y95_l,sp.flipud(y95_u))
    
    y99_l = scale(ls,-sigma99[:,0]+Hmean)
    y99_u = scale(ls,sigma99[:,1]+Hmean)
    y99 = sp.append(y99_l,sp.flipud(y99_u))
    
    p68=Polygon(zip(x,y68),alpha=0.3,lw=0)
    p95=Polygon(zip(x,y95),alpha=0.2,lw=0)
    p99=Polygon(zip(x,y99),alpha=0.1,lw=0)
    
#    plt.figure(1)
#    plt.hold(True)

#    plt.xlabel('$r_{max} [Mpc/h]$')
#    plt.axis([64, 256, 0.90, 1.10])    


    Hmean = scale(ls,Hmean)
    
    plt.gca().add_artist(p68)
    plt.gca().add_artist(p95)
    plt.gca().add_artist(p99)
    plt.plot(ls,Hmean)

    
    return(0)

def solve_qg(parms, q0):

    global cnt
    # Set parameters
    dt = parms.dt
    Nx = parms.Nx
    Ny = parms.Ny

    # initialize fields
    Nt = int(parms.tf/parms.dt)
    #energy = np.zeros(Nt)
    #enstr  = np.zeros(Nt)
    #mass   = np.zeros(Nt)

    # Euler step
    t,ii = 0., 0
    NLnm = flux_qg(q0, parms)
    q  = q0 + dt*NLnm;

    # AB2 step
    t,ii = parms.dt, 1
    NLn = flux_qg(q, parms)
    q   = q + 0.5*dt*(3*NLn - NLnm)

    qe = np.zeros((2*Ny,Nx,2), dtype=float)
    cnt = 1
    for ii in range(3,Nt+1):

        # AB3 step
        t = (ii-1)*parms.dt
        NL = flux_qg(q, parms)
        q  = q + dt/12*(23*NL - 16*NLn + 5*NLnm).real

        # Exponential Filter
        for jj in range(2):
            qe[:,:,jj] = np.vstack((q[:,:,jj],-np.flipud(q[:,:,jj])))
            qe[:,:,jj] = (ifftn(parms.sfilt*fftn(qe[:,:,jj]))).real
            q[:,:,jj]  = qe[0:Ny,:,jj]

        # Reset fluxes
        NLnm = NLn
        NLn  = NL

        if (ii-0)%parms.npt==0:

            # make title
            name = "PV at t = %5.2f" % (t/(3600.0*24.0))

            # Plot PV (or streamfunction)
            plot_q_qhat(q, t)
            #plt.draw()

            cnt += 1

    return q

def elevDeriv(val_mat, direction_mat, inc):

	import scipy;
	import scipy.linalg;
	import math;

	cell_angles  = scipy.flipud(scipy.array([[3 * math.pi / 4, math.pi / 2, math.pi / 4], [math.pi, scipy.nan, 0], [-3 * math.pi / 4, -1 * math.pi / 2, -1 * math.pi / 4]]));
	cell_incs    = scipy.array([[(inc**2 + inc**2)**0.5, inc, (inc**2 + inc**2)**0.5], [inc, scipy.nan, inc], [(inc**2 + inc**2)**0.5, inc, (inc**2 + inc**2)**0.5]]);
	angle        = direction_mat[1,1];

	cell_cosines_f = scipy.cos(angle - cell_angles);
	cell_cosines_b = scipy.cos(angle - cell_angles);
	cell_sines_f   = scipy.sin(angle - cell_angles);
	cell_sines_b   = scipy.sin(angle - cell_angles);

	cell_cosines_f[cell_cosines_f < 0.00001] = scipy.nan;
	cell_cosines_f = cell_cosines_f**2;
	cell_cosines_f = cell_cosines_f / sum(cell_cosines_f[~scipy.isnan(cell_cosines_f)]);
	cell_cosines_b[cell_cosines_b > -0.00001] = scipy.nan;
	cell_cosines_b = cell_cosines_b**2;
	cell_cosines_b = cell_cosines_b / sum(cell_cosines_b[~scipy.isnan(cell_cosines_b)]);
	cell_sines_f[cell_sines_f < 0.00001] = scipy.nan;
	cell_sines_f   = cell_sines_f**2;
	cell_sines_f   = cell_sines_f / sum(cell_sines_f[~scipy.isnan(cell_sines_f)]);
	cell_sines_b[cell_sines_b > -0.00001] = scipy.nan;
	cell_sines_b   = cell_sines_b**2;
	cell_sines_b   = cell_sines_b / sum(cell_sines_b[~scipy.isnan(cell_sines_b)]);

	temp  = val_mat * cell_cosines_f;
	h_x_f = sum(temp[~scipy.isnan(temp)]);

	temp  = val_mat * cell_cosines_b;
	h_x_b = sum(temp[~scipy.isnan(temp)]);

	temp  = val_mat * cell_sines_f;
	h_y_f = sum(temp[~scipy.isnan(temp)]);

	temp  = val_mat * cell_sines_b;
	h_y_b = sum(temp[~scipy.isnan(temp)]);

	h_x = scipy.array([h_x_b, val_mat[1,1], h_x_f]);
	h_y = scipy.array([h_y_b, val_mat[1,1], h_y_f]);

	xs = scipy.array([-1 * int(inc), 0, int(inc)]);
	A  = scipy.vstack([xs, scipy.ones(len(xs))]).T;

#	print("A: " + str(A));
#	print("h_x_b: " + str(h_x_b));
#	print("val_mat[1,1]: " + str(val_mat[1,1]));
#	print("h_x_f: " + str(h_x_f));

	dh_dx, intercept = scipy.linalg.lstsq(A, h_x)[0];
	dh_dy, intercept = scipy.linalg.lstsq(A, h_y)[0];

	return dh_dx, dh_dy;

def flux_qg(q, parms):

    Nx = parms.Nx
    Ny = parms.Ny

    psi = np.zeros((2*Ny,Nx,2),dtype=float)
    u = np.zeros((2*Ny,Nx,2),dtype=float)
    v = np.zeros((2*Ny,Nx,2),dtype=float)
    qe  = np.zeros((2*Ny,Nx,2),dtype=float)
    flux = np.zeros((Ny,Nx,2),dtype=float)
    q_x = np.zeros((2*Ny,Nx,2),dtype=float)
    q_y = np.zeros((2*Ny,Nx,2),dtype=float)
    qe_hat = np.zeros((2*Ny,Nx,2),dtype=complex)
    psie_hat = np.zeros((2*Ny,Nx,2),dtype=complex)

    # - (u + U) q_x - (q_y + Q_y) v
    for ii in range(2):
        # Extend and take FFT
        qe[:,:,ii] = np.vstack((q[:,:,ii],-np.flipud(q[:,:,ii])))
        qe_hat[:,:,ii] = fftn(qe[:,:,ii])

        # Compute gradient of PV
        q_x[:,:,ii] = (ifftn( parms.ikx*qe_hat[:,:,ii])).real
        q_y[:,:,ii] = (ifftn( parms.iky*qe_hat[:,:,ii])).real

    # Compute streamfunction
    psie_hat[:,:,0] = parms.W1[:,:,0]*qe_hat[:,:,0] + parms.W1[:,:,1]*qe_hat[:,:,1]
    psie_hat[:,:,1] = parms.W2[:,:,0]*qe_hat[:,:,0] + parms.W2[:,:,1]*qe_hat[:,:,1]

    for ii in range(2):
        psi[:,:,ii] = (ifftn(psie_hat[:,:,ii])).real

        # Compute physical velocities
        u[:,:,ii] = (ifftn(-parms.iky*psie_hat[:,:,ii])).real
        v[:,:,ii] = (ifftn( parms.ikx*psie_hat[:,:,ii])).real

        # Restrict to physical domain
        #q_x[:,:,ii] = q_x[0:parms.Ny,:,ii]
        #q_y[:,:,ii] = q_y[0:parms.Ny,:,ii]
        #u[:,:,ii]   = u[0:parms.Ny,:,ii]
        #v[:,:,ii]   = v[0:parms.Ny,:,ii]
        #psi[:,:,ii] = psi[0:parms.Ny,:,ii]

        # Compute flux
        flux[:,:,ii] = - (u[0:Ny,:,ii] + parms.U[ii])*q_x[0:Ny,:,ii]- (q_y[0:Ny,:,ii] + parms.Q_y[ii])*v[0:Ny,:,ii]

    # FJP: energy should include potential energy
    #energy = 0.5*np.mean(u**2 + v**2) + np.mean(parms.F*psi**2)
    #enstr  = np.mean(q**2)
    #mass   = np.mean(psi)

    return flux

def coefft(v):
  """Computes Chebyshev coefficients of polynomial represented by collocation
  values in v.  This method is scale invariant; the coefficients of the
  expansion on [a,b] depend on function values in the same way as for [-1,1].
  """
  N = len(v)-1
  if N==0:
    return v
  vrv = sp.concatenate((v,sp.flipud(v[1:-1])))
  U = sp.fft(vrv) / N
  a = sp.ones((N+1,))
  a[[0,-1]] = 0.5
  return a * U[0:N+1]

def filtfilt(b,a,x):
    """
    Filter with given parameters forward and in reverse to eliminate
    phase shifts.
    In addition, initial state is calculated with lfilter_zi and 
    mirror images of the sample are added at end and beginning to
    remove edge effects.
    Must be a one-dimensional array only.
    """
    #For now only accepting 1d arrays
    ntaps=max(len(a),len(b))
    edge=ntaps*3

    if x.ndim != 1:
        raise ValueError, "Filtfilt is only accepting 1 dimension arrays."

    #x must be bigger than edge
    if x.size < edge:
        raise ValueError, "Input vector needs to be bigger than 3 * max(len(a),len(b)."

    if len(a) < ntaps:
        a=scipy.r_[a,scipy.zeros(len(b)-len(a))]

    if len(b) < ntaps:
        b=scipy.r_[b,scipy.zeros(len(a)-len(b))]

    zi=_lfilter_zi(b,a)

    #Grow the signal to have edges for stabilizing 
    #the filter with inverted replicas of the signal
    s=scipy.r_[2*x[0]-x[edge:1:-1],x,2*x[-1]-x[-1:-edge:-1]]
    #in the case of one go we only need one of the extrems 
    # both are needed for filtfilt

    (y,zf)=scipy.signal.lfilter(b,a,s,-1,zi*s[0])

    (y,zf)=scipy.signal.lfilter(b,a,scipy.flipud(y),-1,zi*y[-1])

    return scipy.flipud(y[edge-1:-edge+1])

def solve_qg(parms, q0):

    global cnt
    # Set parameters
    dt = parms.dt
    Nx = parms.Nx
    Ny = parms.Ny

    # initialize fields
    Nt = int(parms.tf/parms.dt)
    energy = np.zeros(Nt)
    enstr  = np.zeros(Nt)
    mass   = np.zeros(Nt)

    # Euler step
    t,ii = 0., 0
    NLnm, energy[0], enstr[0], mass[0] = flux_qg(q0, parms)
    q  = q0 + dt*NLnm;

    # AB2 step
    t,ii = parms.dt, 1
    NLn, energy[1], enstr[1], mass[1] = flux_qg(q, parms)
    q   = q + 0.5*dt*(3*NLn - NLnm)

    kx = fftshift((parms.ikx/parms.ikx[0,1]).real)
    ky = fftshift((parms.iky/parms.iky[1,0]).real)

    cnt = 1
    for ii in range(3,Nt+1):

        # AB3 step
        t = (ii-1)*parms.dt
        NL, energy[ii-1], enstr[ii-1], mass[ii-1] = flux_qg(q, parms)
        q  = q + dt/12*(23*NL - 16*NLn + 5*NLnm).real

        # Exponential Filter
        qe = np.vstack((q,-np.flipud(q)))
        qe = (ifftn(parms.sfilt*fftn(qe))).real
        q  = qe[0:Ny,:]

        # Reset fluxes
        NLnm = NLn
        NLn  = NL

        if (ii-0)%parms.npt==0:

            t = ii*dt
            plot_q_qhat(q, t)
            cnt += 1

    return q, energy, enstr, mass

    def LoadImage(self):
        img = imread(self.photo_data['path']+os.path.sep+self.photo_data['Name'])

        if self.flipud_needed:
            self.axes.imshow(flipud(img), interpolation='nearest',animated=True)
        else:
            self.axes.imshow(img, interpolation='nearest',animated=True)

        self.canvas.draw()
        #store the background boundary box
        self.fig_bg = self.canvas.copy_from_bbox(self.axes.bbox)
        #Save limits
        self.xlims = self.axes.get_xlim()
        self.ylims = self.axes.get_ylim()

 def _add_coi(self, color, data_present=None, fill=False):
     n = len(self.series)
     coi_whole = self.coi * self.dt * sp.hstack((sp.arange((n + 1) / 2), 
                         sp.flipud(sp.arange(n / 2))))
     coi_list = [coi_whole]
     baseline = sp.ones(n) * self.period[-1]
     if data_present is not None:
         for i in range(2, len(data_present) - 1):
             if data_present[i - 1] and (not data_present[i]):
                 coi_list.append(circ_shift(coi_whole, i))
                 baseline[i] = 0
             elif not data_present[i]:
                 baseline[i] = 0
             elif (not data_present[i - 1]) and data_present[i]:
                 coi_list.append(circ_shift(coi_whole, i))
     coi_list.append(baseline)
     coi_line = sp.array(coi_list).min(axis=0)
     coi_line[coi_line == 0] = 1e-4
     x = sp.hstack((self.time, sp.flipud(self.time)))
     y = sp.log2(sp.hstack((coi_line, sp.ones(n) * self.period[-1])))
     if fill:
         plt.fill(x, y, color='black', alpha=0.3)
     else:
         plt.plot(self.time, sp.log2(coi_line), color=color, linestyle=':')

def plot_down_saw_spec_correct():
    plt.close('all')
    rate, in_sig = wavfile.read('saw.wav')
    old_rate = 44100
    new_rate = 22050
    in_sig = sp.float32(in_sig)
    fin = anfft.fft(in_sig)
    nsiz = sp.floor(in_sig.size*new_rate/old_rate)
    nsizh = sp.floor(nsiz/2)
    fout = sp.zeros(nsiz)
    fout = fout + 0j
    fout[0:nsizh] = fin[0:nsizh]
    fout[nsiz-nsizh+1:] = sp.conj(sp.flipud(fout[1:nsizh]))
    f = sp.absolute(fout)
    plt.plot(f[0:f.shape[0]/2])
    plt.savefig('sawdownspec.pdf')