NONSTATIONARITY

The edge effects when moving a window around on a wall of data can be confusing. To illustrate Fortran (Ratfor) technique I wrote a routine nonstat() and indented it so you can see that the 1-axis is handled identically to the 2-axis. First look at the looping instructions which for each axis are one line of initialization and test at the top of the loop, and a line for incrementing at the bottom. The nonstat() routine moves around boxes of any specified size, jumping along the t-axis and x-axis in steps of any specified size.

# slide a window around on a wall of data.  Count times each data point used.
#
subroutine nonstat( n1,n2, w1,w2, k1,k2, count)
integer	n1,n2		# size of data wall
integer	w1,w2		# size of window
integer	k1,k2		# increments for jumping along the wall
integer	s1,s2, e1,e2	# starting and ending points of window on wall
integer	j1,j2		# output math size of array of windows
integer	i1,i2
real	count(n1,n2)	# mathematically it is typically smaller
call zero( n1*n2, count)
j2=1; s2=1; e2=s2+w2-1; while( e2<=n2) {
j1=1; s1=1; e1=s1+w1-1; while( e1<=n1) {
	do i1=1,w1 {
	do i2=1,w2 {
		count(s1+i1-1, s2+i2-1) =  count(s1+i1-1, s2+i2-1) + 1.
		}}
	  j1=j1+1; s1=s1+k1; e1=s1+w1-1
	} j1=j1-1
	  j2=j2+1; s2=s2+k2; e2=s2+w2-1
	} j2=j2-1
return; end

Typically we'll be filtering in each window. The filter output is generally a residual that will be minimized. We won't be running the filter off the edges of the window because that would cause truncation effects preventing us from having small window sizes. So in a given data window the number of residual values is less than the number of data values. The residual has some natural alignment with the data, and we have the problem of splicing together the residuals from each window into a global residual. To do this, we have a space of global residuals on a mesh the size of the data mesh on which we keep a count of times a window residual has added in to the global residual. When everything is finished, the effect of overlapping windows can be compensated by dividing by the count of residuals at each mesh point.