exponential decay script (relaxation response)


A simple linear relaxation response leads to an exponentially decaying impulse response. The response to any input can be found by convolution with this impulse response. The following script will apply this convolution with an exponential using a supplied time-constant parameter: tau.

Select text below and use copy / paste to get the text and save to a file. There must be NO space before the first line. ie “#!” are the first characters in the file.

Save as exp-decay.awk , or change messages to suit.

On Linux/Unix-like systems it will need to be given executable permissions.

#!/bin/awk -f

# convolution integral with exp decay kernel eg Laplace 1/s solution
# usage : ./exp-decay.awk filename <tau=12> <neutral=0>
# optional neutral value is subtracted from input data before integration
# use  OFMT="%6.4f"
# tau is time-const of decay fn,  3 tau window is -95% level; 4=98% ; 5=99.3%
# spin-up:  dataset shortened by window length at begining
# ensure data is equally spaced and continuous !!
# nov 2012  OFMT="%10.8f"

BEGIN { OFMT="%10.8f"
# ARGV[1]=filename; argv[0] is script name, hence ARGC>=1

  if  ( ARGC <2 ) {print " usage : ./exp-decay.awk filename <tau=12> <neutral=0>" ; err=1; exit err; }
  if  ( ARGC >3 ) {neutral=ARGV[3];ARGV[3]=""} else {neutral=0};
  if  ( ARGC >2 ) {tau=ARGV[2];ARGV[2]=""} else {tau=12};

  if (neutral != 0) { print "# subtracting neutral point  const of = ",neutral}

kw=4*tau;  # exp(-x) is approx %5 at 3 time-const, prefer 98% accuracy

# calculate normalised kernel coeffs. kernel is time reversed as required
  for (tot_wt=j=0;j<=kw;j++) {tot_wt+=kwt[-j]=exp(-j/tau) };
  for (j=0;j<=kw;j++) {kwt[-j]/=tot_wt};

# strip off last .xxx part of file name

  for (i=1;i<length(fn);i++){   # if multiple dots in name, build back up except last part
    else basename=basename"."fn[i];

  print "# integrating "ARGV[1]" with exponential decay time const of ",tau
  print "# "ARGV[1]" with exponential decay time const of ",tau >out_file;

($0 !~ /^#/)&&($0 != ""){

  if (ln>kw)
    for (j=-kw;j<=0;j++) {sum+=(ydata[ ln+j ])*kwt[j]}
    print NR,xdata[ln],sum
    print xdata[ln],sum >> out_file;  # spin-up: full integral relates to end date of kernel

  if (err) exit;
  print "#exp decay kernel width = "kw",done"
  if ( neutral !=0) print "#neutral value = "neutral;
  print "#output file = "out_file
#  for (j=-kw;j<=0;j++) print kwt[j];