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

  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"-exp_"tau".dat";

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

($0 !~ /^#/)&&($0 != ""){
  xdata[++ln]=$1;
  ydata[ln]=$2-neutral;

  if (ln>kw)
  {
    sum=0
    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
  }
}

END {
  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];
}


Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s