Nonstationarity: patching

Designing a separate filter for each patch

Recall the prediction-error filter subroutine find_pef() . Given a data plane, this subroutine finds a filter that tends to whiten the spectrum of that data plane. The output is white residual. Now suppose we have a data plane where the dip spectrum is changing from place to place. Here it is natural to apply subroutine find_pef() in local patches. This is done by subroutine find_lopef(). The output of this subroutine is an array of helix-type filters, which can be used, for example, in a local convolution operator loconvol user/gee/lopef.c

    patch_init(dim, npatch, nwall, nwind);
    for (ip=0; ip < np; ip++) {
	bb = aa+ip;

	patch_lop(false, false, n, nw, wall, windata);
	if (NULL != mask) {
	    patch_lop(false, false, n, nw, mask, winmask);
	    for (iw=0; iw < nw; iw++) {
		known[iw] = (winmask[iw] != 0.);
	    }
	    find_mask(nw, known, bb);
	}
	for (mis=iw=0; iw < nw; iw++) {
	    if (!bb->mis[iw]) mis++;
	}
	if (mis > nh) { /* enough equations */
	    find_pef(nw, windata, bb, nh);
	} else if (ip > 1) { /* use last PEF */
	    for (ih=0; ih < nh; ih++) {
		bb->flt[ih] = (bb-1)->flt[ih];
	    }
	}
	patch_close();
    }

We notice that when a patch has fewer regression equations than the filter has coefficients, then the filter is taken to be that of the previous patch.

user/gee/loconvol.c

void loconvol_init(sf_filter aa_in)
/*< initialize with the first filter >*/
{
    aa = aa_in;
}

void loconvol_lop(bool adj, bool add, int nx, int ny, 
		  float *xx, float *yy)
/*< convolve >*/
{
    sf_helicon_init(aa);
    aa++;

    sf_helicon_lop(adj, add, nx, ny, xx, yy);
}

Nonstationarity: patching