next up previous [pdf]

Next: References Up: MULTISCALE, SELF-SIMILAR FITTING Previous: Scale-invariance introduces more fitting

Coding the multiscale filter operator

Equation (3) shows an example where the first output signal is the ordinary one and the second output signal used a filter interlaced with zeros. We prepare subroutines that allow for more output signals, each with its own filter interlace parameter given in the table jump[ns]. Each entry in the jump table corresponds to a particular scaling of a filter axis. The number of output signals is ns and the number of zeros interlaced between filter points for the j-th signal is jump[j]-1.

The multiscale helix filter is defined in module mshelix [*], analogous to the single-scale module helix [*]. A new function onescale extracts our usual helix filter of one particular scale from the multiscale filter.

#ifndef _mshelix_h

typedef struct mshelixfilter {
    int nh, ns;
    float* flt;
    int** lag;
    bool** mis;
    sf_filter one;
} *msfilter;


void onescale(int i, msfilter aa) 
/*< select one scale from multiple scales >*/
    aa->one->lag = aa->lag[i];
    if (NULL != aa->mis) aa->one->mis = aa->mis[i];

We create a multscale helix with module createmshelix [*]. An expanded scale helix filter is like an ordinary helix filter except that the lags are scaled according to a jump.

msfilter createmshelix(int ndim    /* number of dimensions */, 
		       int* nd     /* data size [ndim] */, 
		       int* center /* filter center [ndim] */, 
		       int* gap    /* filter gap [ndim] */, 
		       int ns      /* number of scales */, 
		       int *jump   /* filter scaling [ns] */, 
		       int* na     /* filter size [ndim] */)
/*< allocate and output a multiscale helix filter >*/
    msfilter msaa;
    sf_filter aa;
    int is, ih, nh, id, n123, nb[SF_MAX_DIM];

    aa = createhelix(ndim, nd, center, gap, na);
    nh = aa->nh;

    msaa = msallocate(nh, ns);
    for (is=0; is < ns; is++) {
	for (ih=0; ih < nh; ih++) /* expanded scale lags */
	    msaa->lag[is][ih] = aa->lag[ih]*jump[is]; 

    for (id=0; id < ndim; id++) n123 *= nd[id];    
    msaa->mis = sf_boolalloc2(n123,ns);

    aa = msaa->one;
    for (is=0; is < ns; is++) { /* for all scales */
	onescale( is, msaa); /* extract a filter */  
	aa->mis = NULL;
	for (id=0; id < ndim; id++) nb[id] = na[id]*jump[is];
	bound(ndim, nd, nd, nb, aa); /* set up its boundaries */
	for (id=0; id < n123; id++) 
	    msaa->mis[is][id] = aa->mis[id];  /* save them */
	free (aa->mis);
    return msaa;

First we examine code for estimating a prediction-error filter that is applicable at many scales. We simply invoke the usual filter operator hconest [*] for each scale.

    for (is=0; is < msaa->ns; is++) {

The multiscale prediction-error filter finding subroutine is nearly identical to the usual subroutine find_pef() [*]. (That routine cleverly ignores missing data while estimating a PEF.) To easily extend pef to multiscale filters we replace its call to the ordinary helix filter module hconest [*] by a call to mshconest.

void find_pef(int nd      /* data size */, 
	      float* dd   /* data */, 
	      msfilter aa /* estimated filter */, 
	      int niter   /* number of iterations */) 
/*< estimate PEF >*/
    float *ee;
    int is, id, ns;

    ns = aa->ns;
    ee = sf_floatalloc(nd*ns);
    for (is=0; is < ns; is++) {
	for (id=0; id < nd; id++) {
	    ee[id+is*nd] = dd[id];

    mshconest_init( dd, aa);
    sf_solver(mshconest_lop, sf_cgstep, aa->nh, nd*ns, aa->flt, 
	      ee, niter, "x0", aa->flt, "end");

The purpose of pack(dd,.true.) is to produce the one-dimensional array expected by our solver routines.

Similar code applies to the operator in (4) which is needed for missing data applications. This is like mshconest [*] except the adjoint is not the filter but the input.

    for (is=0; is < aa->ns; is++) {
The multiscale missing-data module msmis2 is just like the usual missing-data module mis2 [*] except that the filtering is done with the multiscale filter mshelicon.

void msmis2(int niter         /* number of iterations */, 
	    int nx            /* data size */, 
	    int ns            /* number of scales */, 
	    float *xx         /* data */, 
	    msfilter aa       /* filter */, 
	    const bool *known /* mask for known data */)
/*< interpolate >*/ 
    int ix, nxs;
    float *dd;

    nxs = ns*nx;

    dd = sf_floatalloc(nxs);
    for (ix=0; ix < nxs; ix++) {
    sf_solver (mshelicon_lop, sf_cgstep, nx, nxs, xx, dd, niter, 
	       "known", known, "x0", xx, "end");

next up previous [pdf]

Next: References Up: MULTISCALE, SELF-SIMILAR FITTING Previous: Scale-invariance introduces more fitting