From e14342e827e2b42e43c006df90c7ad99e5124b3c Mon Sep 17 00:00:00 2001 From: wukong Date: Sun, 4 Nov 2018 10:59:44 -0800 Subject: added gaussian.awk based on hamming script; added pwr(x,p) function to sterling_approx; minor clean up all around; --- fib.awk | 4 ++-- gaussian.awk | 26 ++++++++++++++++++++++++++ hamming.awk | 33 +++++++++++++++++++++------------ hamming.sh | 6 +++--- lin_reg.awk | 21 +++++++++++---------- lin_reg1.awk | 23 +++++++++++++---------- lin_reg2.awk | 26 +++++++++++++------------- quad_reg.awk | 19 ++++++++++--------- sterling_approx.awk | 7 ++++++- 9 files changed, 105 insertions(+), 60 deletions(-) create mode 100644 gaussian.awk 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 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 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) } -- cgit v1.2.3