def set_data(self, x, y, A):
     x = asarray(x).astype(Float32)
     y = asarray(y).astype(Float32)
     A = asarray(A)
     if len(x.shape) != 1 or len(y.shape) != 1\
        or A.shape[0:2] != (y.shape[0], x.shape[0]):
         raise TypeError("Axes don't match array shape")
     if len(A.shape) not in [2, 3]:
         raise TypeError("Can only plot 2D or 3D data")
     if len(A.shape) == 3 and A.shape[2] not in [1, 3, 4]:
         raise TypeError("3D arrays must have three (RGB) or four (RGBA) color components")
     if len(A.shape) == 3 and A.shape[2] == 1:
          A.shape = A.shape[0:2]
     if len(A.shape) == 2:
         if typecode(A) != UInt8:
             A = (self.cmap(self.norm(A))*255).astype(UInt8)
         else:
             A = repeat(A[:,:,NewAxis], 4, 2)
             A[:,:,3] = 255
     else:
         if typecode(A) != UInt8:
             A = (255*A).astype(UInt8)
         if A.shape[2] == 3:
             B = zeros(tuple(list(A.shape[0:2]) + [4]), UInt8)
             B[:,:,0:3] = A
             B[:,:,3] = 255
             A = B
     self._A = A
     self._Ax = x
     self._Ay = y
     self._imcache = None

def epoch2num(e):
    """
    convert an epoch or sequence of epochs to the new date format,
    days since 0001
    """
    spd = 24.*3600.
    return 719163 + asarray(e)/spd

def intwave(wavelet, precision=8):
    """
    intwave(wavelet, precision=8) -> [int_psi, x]
        - for orthogonal wavelets

    intwave(wavelet, precision=8) -> [int_psi_d, int_psi_r, x]
        - for other wavelets

    intwave((function_approx, x), precision=8) -> [int_function, x]
        - for (function approx., x grid) pair

    Integrate *psi* wavelet function from -Inf to x using the rectangle
    integration method.

    wavelet         - Wavelet to integrate (Wavelet object, wavelet name string
                      or (wavelet function approx., x grid) pair)

    precision = 8   - Precision that will be used for wavelet function
                      approximation computed with the wavefun(level=precision)
                      Wavelet's method.

    (function_approx, x) - Function to integrate on the x grid. Used instead
                           of Wavelet object to allow custom wavelet functions.
    """

    if isinstance(wavelet, tuple):
        psi, x = asarray(wavelet[0]), asarray(wavelet[1])
        step = x[1] - x[0]
        return integrate(psi, step), x

    else:
        if not isinstance(wavelet, WAVELET_CLASSES):
            wavelet = wavelet_for_name(wavelet)

        functions_approximations = wavelet.wavefun(precision)
        if len(functions_approximations) == 2:      # continuous wavelet
            psi, x = functions_approximations
            step = x[1] - x[0]
            return integrate(psi, step), x
        elif len(functions_approximations) == 3:    # orthogonal wavelet
            phi, psi, x = functions_approximations
            step = x[1] - x[0]
            return integrate(psi, step), x
        else:                                       # biorthogonal wavelet
            phi_d, psi_d, phi_r, psi_r, x = functions_approximations
            step = x[1] - x[0]
            return integrate(psi_d, step), integrate(psi_r, step), x

def centfrq(wavelet, precision=8):
    """
    centfrq(wavelet, precision=8) -> float
        - for orthogonal wavelets

    centfrq((function_approx, x), precision=8) -> float
        - for (function approx., x grid) pair

    Computes the central frequency of the *psi* wavelet function.

    wavelet         - Wavelet (Wavelet object, wavelet name string
                      or (wavelet function approx., x grid) pair)
    precision = 8   - Precision that will be used for wavelet function
                      approximation computed with the wavefun(level=precision)
                      Wavelet's method.

    (function_approx, xgrid) - Function defined on xgrid. Used instead
                      of Wavelet object to allow custom wavelet functions.
    """

    if isinstance(wavelet, tuple):
        psi, x = asarray(wavelet[0]), asarray(wavelet[1])
    else:
        if not isinstance(wavelet, WAVELET_CLASSES):
            wavelet = wavelet_for_name(wavelet)
        functions_approximations = wavelet.wavefun(precision)

        if len(functions_approximations) == 2:
            psi, x = functions_approximations
        else:
            # (psi, x)   for (phi, psi, x)
            # (psi_d, x) for (phi_d, psi_d, phi_r, psi_r, x)
            psi, x = functions_approximations[1], functions_approximations[-1]

    domain = float(x[-1] - x[0])
    assert domain > 0

    index = argmax(abs(fft(psi)[1:])) + 2
    if index > len(psi) / 2:
        index = len(psi) - index + 2

    return 1.0 / (domain / (index - 1))

def orthfilt(scaling_filter):
    assert len(scaling_filter) % 2 == 0

    scaling_filter = asarray(scaling_filter, dtype=float64)

    rec_lo = sqrt(2) * scaling_filter / sum(scaling_filter)
    dec_lo = rec_lo[::-1]

    rec_hi = qmf(rec_lo)
    dec_hi = rec_hi[::-1]

    return (dec_lo, dec_hi, rec_lo, rec_hi)

def specgram(x, NFFT=256, Fs=2, detrend=detrend_none,
             window=window_hanning, noverlap=128):
    """
    Compute a spectrogram of data in x.  Data are split into NFFT
    length segements and the PSD of each section is computed.  The
    windowing function window is applied to each segment, and the
    amount of overlap of each segment is specified with noverlap

    See pdf for more info.

    The returned times are the midpoints of the intervals over which
    the ffts are calculated
    """
    x = asarray(x)
    assert(NFFT>noverlap)
    if log(NFFT)/log(2) != int(log(NFFT)/log(2)):
       raise ValueError, 'NFFT must be a power of 2'

    # zero pad x up to NFFT if it is shorter than NFFT
    if len(x)<NFFT:
        n = len(x)
        x = resize(x, (NFFT,))
        x[n:] = 0
    

    # for real x, ignore the negative frequencies
    if typecode(x)==Complex: numFreqs=NFFT
    else: numFreqs = NFFT//2+1
        
    windowVals = window(ones((NFFT,),typecode(x)))
    step = NFFT-noverlap
    ind = arange(0,len(x)-NFFT+1,step)
    n = len(ind)
    Pxx = zeros((numFreqs,n), Float)
    # do the ffts of the slices

    for i in range(n):
        thisX = x[ind[i]:ind[i]+NFFT]
        thisX = windowVals*detrend(thisX)
        fx = absolute(fft(thisX))**2
        # Scale the spectrum by the norm of the window to compensate for
        # windowing loss; see Bendat & Piersol Sec 11.5.2
        Pxx[:,i] = divide(fx[:numFreqs], norm(windowVals)**2)
    t = 1/Fs*(ind+NFFT/2)
    freqs = Fs/NFFT*arange(numFreqs)

    return Pxx, freqs, t