frequency axis is in units of “per year” not “year” as labelled.
Equivalent high-pass filter:
Discussion: https://climategrog.wordpress.com/?attachment_id=659
Source Code: ( to get code, use browser copy / paste function within code block below).
#!/bin/awk -f
# usage: ./lanczos_LP.awk filename [50%_attn_period=18] [lobes=3]"
# purpose: convolve 2nd column of file with lanczos windowed sinc function
# this provides low-pass filter with a fast transition from pass-band to stop-band , with minimal ripple in both.
# column 1 data should be evenly spaced without gaps.
# optional second parameter is period of 50% attenuation, as a period expressed in data intervals
# optional third parameter is the number of lobes in impulse response.
# zero pass req period * 1.436 eg 17.23mo or 525 days
# default=18 ; eg. 18 monthly periods : 50% attenuation at 15mo, close to zero pass at 12 mo. for 3 lobe.
# lobes=3 ; use 525 for daily ; lobes=5 ; use 461 for 365 data
# console output comments give pass-band and stop band info. this is not sent to output file.
# More lobes gives sharper transition band but more overshoot distortion / ringing.
# default=3; 3 and 4 are good. ### >=5 may produce significant ringing artefacts ###
function sinc(x) {
if (x==0) return 1; else return (sin(x)/x);
}
BEGIN{ OFMT = "%8.6f"
if ( ARGC >3 ) {lobes=ARGV[3];ARGV[3]=""} else {lobes=3};
if ( ARGC >2 ) {period=ARGV[2];ARGV[2]=""} else {period=18};
if ( ARGC <2) { print "### usage: ./lanczos_LP.awk filename [50%_attn_period=18] [lobes=3]"; err=1; exit 1; }
else {len=ARGV[2]; ARGC=2;}
pi= 3.14159265359811668006 ;
twopi=2*pi;
halfk=int(period/2*lobes);
kwid= 2*halfk; # kernel = -halfk....+halfk = kwid+1 pts ie odd number
w_taper=twopi/kwid; # first lanczos lobe -pi to +pi to fill window
w=twopi/period; # sinc passes zero pi, 2*pi, 3*pi ....
# calculate normalised kernel coeffs
for (tot_wt=j=0;j<=halfk;j++) {tot_wt+=gwt[-j]=gwt[j]=sinc(j*w)*sinc(j*w_taper)};
tot_wt=2*tot_wt-gwt[0];
for (j=-halfk;j<=halfk;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,".");
for (i=1;i<length(fn);i++){ # if multiple dots in name, build back up except last part
if(i==1)basename=fn[i]
else basename=basename"."fn[i];
}
out_file=basename"-lanc"lobes"LP_"period".dat";
print "# ",out_file >out_file;
print "# ",out_file
}
ln=-1; # init linenumber counter
}
($0 !~ /^#/)&&($0 != "") { # strip out comments, headers and blank lines
xdata[++ln]=$1;
ydata[ln]=$2;
if (ln>2*halfk)
{
sum=0
for (j=-2*halfk;j<=0;j++) {sum+=ydata[ln+j]*gwt[j+halfk]}
if (out_file=="-")
print xdata[ln-halfk],sum
else {
print NR,xdata[ln-halfk],sum
print xdata[ln-halfk],sum >> out_file;
}
} #endif ln
}
END {
if (err) exit err;
if (FNR<=kwid) {print " ### insufficient data to fill "kwid+1" point filter buffer"; exit 1}
print "# "lobes"-lobe lanczos kernel width = "2*halfk+1",done"
if (lobes==5) printf ("# "lobes"-lobe lanczos LP on %s: zero-pass at %.3f; 99% pass > %.3f \n",ARGV[1],period/1.258,period/.752)
if (lobes==4) printf ("# "lobes"-lobe lanczos LP on %s: zero-pass at %.3f; 99% pass > %.3f \n",ARGV[1],period/1.325,period/.698)
if (lobes==3) printf ("# "lobes"-lobe lanczos LP on %s: zero-pass at %.3f; 99% pass > %.3f \n",ARGV[1],period/1.436,period/.588)
print "# output file = "out_file
# uncomment following to dump kernel (impluse response):
# sum=0; for (j=-halfk;j<=halfk;j++) {sum+= gwt[j]; printf ("%.4f\n", gwt[j])}
# print "# abs sum of coeffs = ",sum;
}
climategrog 