Spatial aliasing and scale invariance

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.

user/gee/mshelix.c

#ifndef _mshelix_h

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

#endif

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.

user/gee/createmshelix.c

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]; 
	
    }
    sf_deallocatehelix(aa);

    n123=1;
    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.

user/gee/mshconest.c

    for (is=0; is < msaa->ns; is++) {
	onescale(is,msaa);
	hconest_lop(adj,true,na,nx,a,y+is*nx);
    }

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.

user/gee/mspef.c

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");
    sf_cgstep_close();

    free(ee);
}

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.

user/gee/mshelicon.c

    for (is=0; is < aa->ns; is++) {
	onescale(is,aa);
	sf_helicon_init(aa->one);
	sf_helicon_lop(adj,true,nx,nx,xx,yy+is*nx);
    }

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.

user/gee/msmis2.c

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++) {
	dd[ix]=0.;
    }
    
    mshelicon_init(aa);
    sf_solver (mshelicon_lop, sf_cgstep, nx, nxs, xx, dd, niter, 
	       "known", known, "x0", xx, "end");
    
    sf_cgstep_close();
    free(dd);
}

Spatial aliasing and scale invariance