added gaussian.awk based on hamming script;

added pwr(x,p) function to sterling_approx;
minor clean up all around;

9 files changed, 105 insertions, 60 deletions

diff --git a/fib.awk b/fib.awk
index 07e8bab..d599caa 100644
--- a/fib.awk
+++ b/fib.awk
@@ -5,6 +5,6 @@
 
 BEGIN {
     n = ARGV[1]
-    printf(OFMT ORS,
-        (1.0/sqrt(5.0))*((1.0 + sqrt(5.0))/2.0)^n  - (1.0/sqrt(5.0))*((1.0 - sqrt(5.0))/2.0)^n)
+    print (1/sqrt(5))*(0.5*(1 + sqrt(5)))^n - (1/sqrt(5))*(0.5*(1 - sqrt(5)))^n
 }
+
diff --git a/gaussian.awk b/gaussian.awk
new file mode 100644
index 0000000..ef81b68
--- /dev/null
+++ b/gaussian.awk
@@ -0,0 +1,26 @@
+#!/usr/bin/awk -f

+

+### gaussian.awk

+# generate a Gaussian window

+# https://en.wikipedia.org/wiki/Window_function

+

+BEGIN {

+

+    N = ARGV[1]

+    M = 0.5*(N - 1)

+

+    # sigma <= 0.5

+    sigma = 0.4

+

+    for (n=0; n<N; n++) {

+        if (N > 1 && M > 0) {

+            e[n] = (n - M)/(sigma*M)

+            w[n] = exp(-0.5*e[n]*e[n])

+            print n, w[n]/M

+        }

+        else {

+            print n, 1.0

+        }

+    }

+

+}

diff --git a/hamming.awk b/hamming.awk
index 95f1cde..fd6d89c 100644
--- a/hamming.awk
+++ b/hamming.awk
@@ -6,25 +6,34 @@
 # 'a0' and 'a1' parameters of the raised cosine.

 

 BEGIN {

+

+    ### R.W. Hamming, "Digital Filters"

+    # H = "0.23 0.54 0.23"

+

     N = ARGV[1]

-    pi = 4*(4*atan2(1,5) - atan2(1,239))

+    # window interval goes from -M to M

+    M = 0.5*(N - 1)

+    pi = 4*atan2(1,1)

 

-    #a0 = (25.0/46.0)

-    #a1 = (21.0/46.0)

+    a0 = (25.0/46.0)

+    a1 = (21.0/46.0)

 

     ### optimal values for equal-ripple

-    a0 = 0.53836

-    a1 = 0.46164

+    #a0 = 0.53836

+    #a1 = 0.46164

 

-    ### R.W. Hamming, "Digital Filters"

-    # H = "0.23 0.54 0.23"

+    ### VonHann

+    #a0 = 0.5

+    #a1 = 0.5

 

-    for (n=0; n<N; n++) {

-        if (N > 1) {

-            w[n] = a0 - a1*cos((2*pi*n)/(N - 1))

-            sum_w += w[n]

+    for (n=-M; n<=M; n++) {

+        if (N > 1 && M > 0) {

+            w[n] = a0 + a1*cos((pi*n)/M)

+            print n, w[n]/M

+        }

+        else {

+            print n, 1.0

         }

-        print 2.0*(w[n])/(N + 1)

     }

 

 }

diff --git a/hamming.sh b/hamming.sh
index fdd6ba6..b466325 100644
--- a/hamming.sh
+++ b/hamming.sh
@@ -2,9 +2,9 @@
 
 ### hamming.sh
 # Convolve input vector with 3-element Hamming window [0.23, 0.54, 0.23]
-# note: bash built-in $(( )) constructs only support integer arithmetic
-# floating point operation will require another language such as awk, perl, or
-# python.
+# note: bash built-in $(( )) construct only supports integer arithmetic
+# floating point operations will require another language such as awk, perl,
+# python, etc.
 
 ### input vector as array
 X=($1)
diff --git a/lin_reg.awk b/lin_reg.awk
index e0d67ca..205c43f 100644
--- a/lin_reg.awk
+++ b/lin_reg.awk
@@ -37,7 +37,7 @@ NF > 0 {
             sum_delta2[y] += delta0[y]*delta1[y]
 
             ### sample variance
-            (count[y] > 1) ? var[y] = sum_delta2[y]/(count[y] - 1) : var[y] = ""
+            #(count[y] > 1) ? var[y] = sum_delta2[y]/(count[y] - 1) : var[y] = ""
 
             # x = row, y = col, trendline: y = A + Bx
             for (x=1; x<=nf_max; x++) {
@@ -46,7 +46,12 @@ NF > 0 {
                 sum_delta_xy[x,y] += delta0[x]*delta1[y]
 
                 # covariance
-                (count[x,y] > 1) ? cov[x,y] = sum_delta_xy[x,y]/(count[x,y] - 1) : cov[x,y] = ""
+                #(count[x,y] > 1) ? cov[x,y] = sum_delta_xy[x,y]/(count[x,y] - 1) : cov[x,y] = ""
+
+
+                # correlation
+                r_den[x,y] = sqrt(sum_delta2[x]*sum_delta2[y])
+                (r_den[x,y]) ? r[x,y] = sum_delta_xy[x,y]/r_den[x,y] : r[x,y] = 1
 
                 ab_den[x,y] = (count[x,y]*sum2[x] - sum[x]*sum[x])
                 if (ab_den[x,y]) {
@@ -70,11 +75,6 @@ NF > 0 {
                 a_err_den[x,y] = count[x,y]*b_err_den[x,y]
                 if (a_err_den[x,y])
                     a_err[x,y] = sqrt(sum2[x]/count[x,y])*b_err[x,y]
-
-                # correlation
-                r_den[x,y] = sqrt(sum_delta2[x]*sum_delta2[y])
-                (r_den[x,y]) ? r[x,y] = sum_delta_xy[x,y]/r_den[x,y] : r[x,y] = 1
-
             }
         }
         else
@@ -86,9 +86,10 @@ END {
     for (y=1; y<=nf_max; y++) {
         for (x=1; x<=nf_max; x++) {
             if (x != y && r[x,y]) {
-                printf(OFMT OFS "(%s)" OFS " = (" OFMT " +/- " OFMT ")(%s)" OFS " + (" OFMT " +/- "OFMT ")%s",
-                    10.0*log(r[x,y]*r[x,y])/log(10.0), header[y], b[x,y],
-                    b_err[x,y], header[x], a[x,y], a_err[x,y], ORS)
+                printf(OFMT OFS "(%s)" OFS " = (" OFMT " +/- " OFMT ")(%s)" OFS " + (" OFMT " +/- " OFMT ")" ORS,
+                    10.0*log(r[x,y]*r[x,y])/log(10.0),
+                    header[y], b[x,y], b_err[x,y],
+                    header[x], a[x,y], a_err[x,y])
             }
         }
     }
diff --git a/lin_reg1.awk b/lin_reg1.awk
index 9d9f257..ae59c77 100644
--- a/lin_reg1.awk
+++ b/lin_reg1.awk
@@ -22,28 +22,30 @@ NF > 0 {
 
     ### iterate over columns
     for (y=1; y<=nf_max; y++) {
+        if ($y == header[n])
+            continue
         if ($y ~ number) {
 
             ### mean
             count[y] += 1
             sum[y] += $y
             sum2[y] += $y*$y
-            mean[y] = sum[y]/count[y]
+            delta0[y] = $y - mean[y]
+            mean[y] += delta0[y]/count[y]
+            delta1[y] = $y - mean[y]
+            sum_delta[y] += delta1[y]
+            sum_delta2[y] += delta0[y]*delta1[y]
 
-            ### difference from the mean
-            delta[y] = $y - mean[y]
-            sum_delta[y] += delta[y]
-            sum_delta2[y] += delta[y]*delta[y]
 
-            # x = row, y = col
+            # x = row, y = col, trendline: y = A + Bx
             for (x=1; x<=nf_max; x++) {
                 count[x,y] += 1
                 sum_xy[x,y] += $x*$y
-                sum_delta_xy[x,y] += delta[x]*delta[y]
+                sum_delta_xy[x,y] += delta0[x]*delta1[y]
 
                 # correlation
                 r_den[x,y] = sqrt(sum_delta2[x]*sum_delta2[y])
-                (r_den[x,y]) ? r[x,y] = sum_delta_xy[x,y]/r_den[x,y] : r[x,y] = 0
+                (r_den[x,y]) ? r[x,y] = sum_delta_xy[x,y]/r_den[x,y] : r[x,y] = 1
 
                 ab_den[x,y] = (count[x,y]*sum2[x] - sum[x]*sum[x])
                 if (ab_den[x,y]) {
@@ -65,9 +67,10 @@ END {
     for (y=1; y<=nf_max; y++) {
         for (x=1; x<=nf_max; x++) {
             if (x != y && r[x,y]) {
-                r2[x,y] = r[x,y]*r[x,y]
                 printf(OFMT OFS "(%s)" OFS " = " OFMT "(%s)" OFS " + " OFMT ORS,
-                    10.0*log(r2[x,y])/log(10.0), header[y], b[x,y], header[x], a[x,y])
+                    10.0*log(r[x,y]*r[x,y])/log(10.0),
+                    header[y], b[x,y],
+                    header[x], a[x,y])
             }
         }
     }
diff --git a/lin_reg2.awk b/lin_reg2.awk
index 8939947..f06d890 100644
--- a/lin_reg2.awk
+++ b/lin_reg2.awk
@@ -22,27 +22,26 @@ NF > 0 {
 
     ### iterate over columns
     for (y=1; y<=nf_max; y++) {
+        if ($y == header[n])
+            continue
         if ($y ~ number) {
 
             ### mean
             count[y] += 1
             sum[y] += $y
             sum2[y] += $y*$y
-            mean[y] = sum[y]/count[y]
-
-            ### difference from the mean
-            delta[y] = $y - mean[y]
-            sum_delta[y] += delta[y]
-            sum_delta2[y] += delta[y]*delta[y]
+            delta0[y] = $y - mean[y]
+            mean[y] += delta0[y]/count[y]
+            delta1[y] = $y - mean[y]
+            sum_delta[y] += delta1[y]
+            sum_delta2[y] += delta0[y]*delta1[y]
 
-            ### sample variance
-            (count[y] > 1) ? var[y] = sum_delta2[y]/(count[y] - 1) : var[y] = ""
 
-            # x = row, y = col
+            # x = row, y = col, trendline: y = A + Bx
             for (x=1; x<=nf_max; x++) {
                 count[x,y] += 1
                 sum_xy[x,y] += $x*$y
-                sum_delta_xy[x,y] += delta[x]*delta[y]
+                sum_delta_xy[x,y] += delta0[x]*delta1[y]
 
                 # correlation
                 r_den[x,y] = sqrt(sum_delta2[x]*sum_delta2[y])
@@ -58,7 +57,7 @@ NF > 0 {
                     b[x,y] = 1
                 }
 
-                ### error estimate
+                # error estimate
                 err_den[x,y] = count[x,y]*(count[x,y] - 2)
                 if (count[x,y] > 2) {
                     err[x,y] = $y - (a[x,y] + b[x,y]*$x)
@@ -82,8 +81,9 @@ END {
         for (x=1; x<=nf_max; x++) {
             if (x != y && r[x,y]) {
                 printf(OFMT OFS "(%s)" OFS " = (" OFMT " +/- " OFMT ")(%s)" OFS " + (" OFMT " +/- " OFMT ")" ORS,
-                    10.0*log(r[x,y]*r[x,y])/log(10), header[y], b[x,y], b_err[x,y], header[x],
-                    a[x,y], a_err[x,y])
+                    10.0*log(r[x,y]*r[x,y])/log(10.0),
+                    header[y], b[x,y], b_err[x,y],
+                    header[x], a[x,y], a_err[x,y])
             }
         }
     }
diff --git a/quad_reg.awk b/quad_reg.awk
index 4267412..ede06a4 100644
--- a/quad_reg.awk
+++ b/quad_reg.awk
@@ -32,18 +32,16 @@ NF > 0 {
             sum2[y] += $y*$y
             sum3[y] += $y*$y*$y
             sum4[y] += $y*$y*$y*$y
-            mean[y] = sum[y]/count[y]
-            mean2[y] = sum2[y]/count[y]
-
-            ### delta, difference from the mean
-            delta[y] = $y - mean[y]
-            delta2[y] = $y*$y - mean2[y]
+            delta0[y] = $y - mean[y]
+            mean[y] += delta0[y]/count[y]
+            delta1[y] = $y - mean[y]
+            delta[y] = delta1[y]
+            delta2[y] = delta0[y]*delta[1]
             sum_delta[y] += delta[y]
-            sum2_delta[y] += delta[y]*delta[y]
             sum_delta2[y] += delta2[y]
 
             ### sample variance
-            (count[y] > 1) ? var[y] = sum_delta2[y]/(count[y] - 1) : var[y] = ""
+            #(count[y] > 1) ? var[y] = sum_delta2[y]/(count[y] - 1) : var[y] = ""
 
             # x = row, y = col, trendline: y = A + Bx + Cx^2
             for (x=1; x<=nf_max; x++) {
@@ -122,7 +120,10 @@ END {
                 printf(OFMT OFS "(%s)" OFS " = (" OFMT ")(%s)^2" OFS " + (" OFMT ")(%s)" OFS " + (" OFMT ")",
                     10.0*log(r[x,y]*r[x,y])/log(10), header[y], c[x,y], header[x], b[x,y], header[x], a[x,y])
                 printf("  [" OFMT "," OFMT "][" OFMT "," OFMT "][" OFMT "," OFMT "]" OFS" [" OFMT "," OFMT "]",
-                    rx0[x,y], 0, rx1[x,y], 0, 0, a[x,y], xv[x,y], yv[x,y])
+                    rx0[x,y], (a[x,y] + b[x,y]*rx0[x,y] + c[x,y]*rx0[x,y]*rx0[x,y]),
+                    rx1[x,y], (a[x,y] + b[x,y]*rx1[x,y] + c[x,y]*rx1[x,y]*rx1[x,y]),
+                    0, a[x,y],
+                    xv[x,y], yv[x,y])
                 printf("  [" OFMT "," OFMT "]", xf[x,y], yf[x,y])
                 printf("  [" OFMT "]" ORS, sqrt((yf[x,y] - yv[x,y])*(yf[x,y] - yv[x,y])))
             }
diff --git a/sterling_approx.awk b/sterling_approx.awk
index ca119ac..62825b3 100644
--- a/sterling_approx.awk
+++ b/sterling_approx.awk
@@ -12,12 +12,17 @@
 ### sterling_approx.awk
 # https://en.wikipedia.org/wiki/Stirling%27s_approximation
 
+function pwr(x, p) {
+    return p ? exp(p*log(x)) : 1
+}
+
 BEGIN {
+    #OFMT = "%.18g"
     ARGV[1] ? n = ARGV[1] : n = 0
     pi = 4*atan2(1,1)
     f = 0
     if (n > 0) {
-        f = sqrt(2*pi*n)*exp(n*log(n*exp(-1)))*(1 + 1/(12*n) + 1/(288*n*n) - 139/(51840*n*n*n) - 571/(2488320*n*n*n*n))
+        f = sqrt(2*pi*n)*pwr(n*exp(-1), n)*(1 + 1/(12*n) + 1/(288*n*n) - 139/(51840*n*n*n) - 571/(2488320*n*n*n*n))
     }
     printf(OFMT ORS, f)
 }