undec-interp-comparison.ipynb

undecimated.py

import numpy as np

from scipy.signal import fftconvolve
from scipy.ndimage import maximum_filter
from dtcwt.backend.backend_numpy.lowlevel import reflect

import dtcwt.backend.backend_numpy as dtcwt_backend

def conv(X, h):
    # Extend X reflecting appropriately
    m = tuple(int(np.fix(x*0.5)) for x in h.shape)
    row_idxs = reflect(np.arange(-m[0],X.shape[0]+m[0]-1), -0.5, X.shape[0]-0.5)
    col_idxs = reflect(np.arange(-m[1],X.shape[1]+m[1]-1), -0.5, X.shape[1]-0.5)
    X = (X[row_idxs, :])[:, col_idxs]
    return fftconvolve(X, h, 'valid')

def undec_transform(X, nlevels=3, biort='near_sym_b', qshift='qshift_d', scales=None):
    t = dtcwt_backend.Transform2d(biort=biort, qshift=qshift)

    sz = 2 + (4<<nlevels)
    im = np.zeros((sz,sz))
    ft = t.forward(im, nlevels=nlevels)

    # A sequence of complex impulse responses
    impulse_responses = []
    #scales = [1j, -1j, 1j, -1, 1, -1]

    for subband_idx in xrange(6):
        # Copy the zero result. The ugliness of this implies that I should probably add
        # a .copy() method to the TransformDomainSignal interface.
        mt = dtcwt_backend.TransformDomainSignal(ft.lowpass, list(np.copy(x) for x in ft.subbands))
        sbs = mt.subbands[-1].shape[:2]
        
        mt.subbands[-1][sbs[0]>>1,sbs[1]>>1,subband_idx] = 1
        real_ir = t.inverse(mt).value
        
        mt.subbands[-1][sbs[0]>>1,sbs[1]>>1,subband_idx] = 1j
        imag_ir = t.inverse(mt).value
        
        if scales is not None:
            s = scales[subband_idx]
        else:
            s = 1

        impulse_responses.append((real_ir + 1j*imag_ir) * s)

    subbands = []
    for idx in xrange(len(impulse_responses)):
        subbands.append(conv(X, impulse_responses[idx].real) + 1j * conv(X, impulse_responses[idx].imag))
    subbands = np.dstack(subbands)

    low = conv(X, np.abs(impulse_responses[0]))

    return low, subbands, impulse_responses