The graph is scaled to unit frequency (eg 1 Hz or one cycle per year …).

It illustrates how a 159 day gaussian will filter daily data. For monthly data 12*159/365=5.23 , so the nearest would be sigma=5mo . However since that is slightly short of the intended value and gaussian leaks a little, since it is never fully zero, six months would be better to suppress a 12mo cycle.

Code ( use copy / paste within code block ).

There must be NO space before the first line. ie “#!” are the first characters in the file.

#!/bin/awk -f # pass input through 3 sigma gaussian filter where sigma, if not given, is 2 data points wide # usage : ./gauss.awk filename <sigma=2> <scale_factor=1> # optional scale_factor simply scales the output # sigma can be compared to the period of the -3dB point of the filter # result is centred, ie not shift. dataset shortened by half window each end # operates on files with two columns of numeric data, or a single column of data # eg date as integer year or floating point decimal # single column takes line number position to be evenly spaced ordinate variable. # ensure that data is equally spaced and continuous !! # nov 2012 OFMT="%10.8f" # 159d gauss similar to 365d r3m # FWHM (-3dB) = 2*sigma*sqrt(2*ln2) # half height freq = sigma * 1.1774 # 3 sigma gives 99.7% accuracy BEGIN { OFMT="%10.8f" # ARGV[1]=filename; argv[0] is script name, hence ARGC>=1 if (ARGV[1]=="") { print " # usage : ./gauss.awk filename <sigma=2> <scale_factor=1>" ; exit 1} pi= 3.14159265359811668006 if ( ARGC >3 ) {scaleby=ARGV[3];ARGV[3]=""} else {scaleby=1}; if ( ARGC >2 ) {sigma=ARGV[2];ARGV[2]=""} else {sigma=2}; print "# filtering "ARGV[1]" with gaussian of sigma= ",sigma root2pi_sigma=sqrt(2*pi)*sigma; two_sig_sqr=2.0*sigma*sigma; gw=3*sigma-1; # gauss is approx zero at 3 sigma, use 3 sig window # calculate normalised gaussian coeffs for (tot_wt=j=0;j<=gw;j++) {tot_wt+=gwt[-j]=gwt[j]=exp(-j*j/two_sig_sqr)/ root2pi_sigma}; tot_wt=2*tot_wt-gwt[0]; tot_wt/=scaleby; for (j=-gw;j<=gw;j++) {gwt[j]/=tot_wt}; # strip off last .xxx part of file name # improve this (doesn't work with paths like ../filename) if (ARGV[1]=="-") { out_file="-"; } else { split(ARGV[1],fn,"."); basename=fn[1] out_file=basename"-gauss"sigma".dat"; print "# ",out_file >out_file; } ln=-1; } ( ($0 !~ /^#/) && ($0 !~ /^[" ",\t]*$/) ) { xdata[++ln]=$1; ydata[ln]=$2; if (ln>2*gw) { gauss=0 for (j=-2*gw;j<=0;j++) {gauss+=ydata[ln+j]*gwt[j+gw]} if (out_file=="-") print xdata[ln-gw],gauss else { print NR,xdata[ln-gw],gauss print xdata[ln-gw],gauss >> out_file; } } } END { if (ARGV[1]=="") { exit 1} print "# gausssian window width = "gw+gw+1",done" print "# output file = "out_file # for (i=0;i<=gw;i++){print gwt[i]} }