diff --git a/.gitignore b/.gitignore
index 204659a67114a94a31408bd0b57dfb453ee124c8..14bc3761fb04d7b063e5d32f2937d54e71c3a6e0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,3 +57,7 @@ dkms.conf
# End of https://www.gitignore.io/api/c
lpsd-exec
+example/test.dat
+example/lpsd_test
+example/*.gnu
+example/*.txt
diff --git a/Makefile b/Makefile
index 87be461bdce515476118923d95ff6188fd7ebe67..ab1e7afe45eeea386e74d7706f04073ad02f6d41 100755
--- a/Makefile
+++ b/Makefile
@@ -21,6 +21,16 @@ install:
.PHONY: clean install
+example-data:
+ @echo "Compiling lpsd_test and creating example file test.dat."
+ make -C example/
+ @echo "Running lpsd with test file..."
+ # ./lpsd-exec --input=example/test.dat
+
+clean-example-data:
+ make -C example/ clean
+
clean:
rm $(OBJECTS) lpsd-exec
+
diff --git a/README.md b/README.md
index a367fcb57badc54e0783c51ba4a2e5730b67f4fd..8cc6a60c1480490bccf4161f3e69ff18682e05d9 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,7 @@
-# LPSD Algorithm
+# LPSD Algorithm
+
+> lpsd - a program to calculate spectral estimates by Michael Troebs and
+> Gerhard Heinzel, (c) 2003, 2004
This repository provides the implementation of the LPSD algorithm for spectrum
and spectral density estimation from long time series by Fourier transforms at
@@ -8,14 +11,20 @@ It is written in C by Michael Tröbs and [Gerhard Heinzel](mailto:gerhard.heinze
## Table of Contents
+<!-- vim-markdown-toc GFM -->
+
* [Install](#install)
- * [Requirements](#requirements)
- * [Compilation](#compilation)
+ * [Requirements](#requirements)
+ * [Compilation](#compilation)
* [Usage](#usage)
+ * [Example](#example)
+ * [Options](#options)
* [Documentation](#documentation)
- * [Theoretical Background](#theoretical-background)
- * [File Overview](#file-overview)
- * [Options](#options)
+ * [Theoretical Background](#theoretical-background)
+ * [File Overview](#file-overview)
+* [FAQ](#faq)
+
+<!-- vim-markdown-toc -->
## Install
@@ -25,11 +34,76 @@ Make sure you have `fftw3` installed.
### Compilation
-Make use of the `Makefile` by typing `make`.
+Make use of the `Makefile` to compile the `lpsd-exec` by typing:
+
+```
+$ make
+```
## Usage
-`lpsd` can be controlled by command line options or interactively.
+`lpsd` can be controlled by command line options or interactively.
+
+```
+$ ./lpsd-exec [OPTION...]
+```
+
+Run `lpsd-exec` with the `help` option for a good overview of lpsd's
+functionality and options:
+
+```
+$ ./lpsd-exec --help
+```
+
+### Example
+
+Create a test dataset with:
+
+```
+$ make example-data
+```
+
+This will create the file `test.dat` in the directory `example` which can be
+processed with `lpsd`:
+
+```
+$ ./lpsd-exec --input=example/test.dat
+```
+
+### Options
+
+The command options `lpsd` understands:
+
+| Short | Long | Description |
+| :---: | :----------------------- | :---------------------------------------------- |
+| `-a` | `--davg=# of des. avgs` | desired number of averages |
+| `-A` | `--colA=# of column ` | number of column to process |
+| `-b` | `--tmin=tmin ` | start time in seconds |
+| `-B` | `--colB=# of column ` | process column B - column A |
+| `-c` | `--param=param ` | parameter string |
+| `-d` | `--usedefs ` | use defaults |
+| `-e` | `--tmax=tmax ` | stop time in seconds |
+| `-f` | `--fsamp=sampl. freq. ` | sampling frequency in Hertz |
+| `-g` | `--gnuplot=gnuplot file` | gnuplot file name |
+| `-h` | `--method=0, 1 ` | method for frequency calculation: 0-LPSD, 1-FFT |
+| `-i` | `--input=input file ` | input file name |
+| `-j` | `--fres=FFT freq. res.` | Frequency resolution for FFT |
+| `-k` | `--sbin=sbin ` | smallest frequency bin |
+| `-l` | `--ovlp=overlap ` | segment overlap in % |
+| `-m` | `--mavg=# of min. avgs ` | minimum number of averages |
+| `-n` | `--nspec=# in spectr.` | number of values in spectrum |
+| `-o` | `--output=output file ` | output file name |
+| `-p` | `--psll=psll ` | peak side lobe level in dB |
+| `-q` | `--quiet ` | Don't produce output on screen |
+| `-r` | `--lr=0,1 ` | linear regression; 1 yes, 0 no |
+| `-s` | `--fmin=fmin ` | start frequency in spectrum |
+| `-t` | `--fmax=fmax ` | stop frequency in spectrum |
+| `-T` | `--time ` | file contains time in s in first column |
+| `-u` | `--gnuterm=gnuterm ` | number of gnuplot terminal |
+| `-w` | `--window=wind. func.` | window function; -2 Kaiser, -1 flat top, 0..30 |
+| `-x` | `--scale=factor ` | scaling factor |
+| `-?` | `--help ` | Give this help list |
+| `-V` | `--version ` | Print program version |
## Documentation
@@ -65,38 +139,6 @@ For a more theoretical understanding read the following publication:
| `StrParser.c` | |
| `tics.c` | |
-### Options
-
-The options `lpsd` understands:
-
-| Option | Short | Value | Description |
-| :------: | :---: | :------------- | ----------------------------------------------- |
-| `davg` | `a` | # of des. avgs | desired number of averages |
-| `colA` | `A` | # of column | number of column to process |
-| `tmin` | `b` | tmin | start time in seconds |
-| `colB` | `B` | # of column | process column B - column A |
-| `param` | `c` | param | parameter string |
-| `usedefs` | `d` | | use defaults |
-| `tmax` | `e` | tmax | stop time in seconds |
-| `fsamp` | `f` | sampl. freq. | sampling frequency in Hertz |
-| `gnuplot` | `g` | gnuplot file | gnuplot file name |
-| `method` | `h` | 0, 1 | method for frequency calculation: 0-LPSD, 1-FFT |
-| `input` | `i` | input file | input file name |
-| `fres` | `j` | FFT freq. res. | Frequency resolution for FFT |
-| `sbin` | `k` | sbin | smallest frequency bin |
-| `ovlp` | `l` | overlap | segment overlap in % |
-| `mavg` | `m` | # of min. avgs | minimum number of averages |
-| `nspec` | `n` | # in spectr. | number of values in spectrum |
-| `output` | `o` | output file | output file name |
-| `psll` | `p` | psll | peak side lobe level in dB |
-| `quiet` | `q` | 0 | Don't produce output on screen |
-| `lr` | `r` | 0,1 | linear regression; 1 yes, 0 no |
-| `fmin` | `s` | fmin | start frequency in spectrum |
-| `fmax` | `t` | fmax | stop frequency in spectrum |
-| `time` | `T` | | tile contains time in s in first column |
-| `gnuterm` | `u` | gnuterm | number of gnuplot terminal |
-| `window` | `w` | wind. func. | window function; -2 Kaiser, -1 flat top, 0..30 |
-
## FAQ
* Hot can I enable `lpsd` configuration via config file?
@@ -125,5 +167,3 @@ Example:
```
$ set LPSDCFN=c:\tools\lpsd.cfg
```
-
-
diff --git a/example/Makefile b/example/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..1095f4d571269ffd4409e9a0aa6026c8833b16d5
--- /dev/null
+++ b/example/Makefile
@@ -0,0 +1,9 @@
+
+all:
+ gcc -O -Wall -W lpsd_test.c random.c -o lpsd_test -lm
+ ./lpsd_test 1000000 > test.dat
+
+.PHONY: clean
+
+clean:
+ rm -f lpsd_test test.dat
diff --git a/example/lpsd_test.c b/example/lpsd_test.c
new file mode 100644
index 0000000000000000000000000000000000000000..0c84a92c1b054905fc1255da2b6243ef4f9012c4
--- /dev/null
+++ b/example/lpsd_test.c
@@ -0,0 +1,407 @@
+/*
+gcc -O -Wall -W lpsd_test.c random.c -o lpsd_test -lm -lgmp && ./lpsd_test 50
+*/
+
+/* code generated by
+LISO analog circuit simulation 2.06b (Aug 28 2017) gerhard.heinzel@aei.mpg.de
+Input file lpsd_test1.fil, Output file lpsd_test1.out, Gnuplot file lpsd_test1.gnu
+ghh@hws393 Mon Mar 5 12:46:40 2018
+*/
+
+/*
+ pole 100 3
+ zero 300 0.7
+ pole 2000
+ rndgenalias 20k
+ freq log 1 10k 1000
+
+ LISO internal parameters:
+ internal MPF precision: 2048 bits (RND_MPF_PREC) [approx. 616 decimal digits]
+ convergence criterion for Matrix exponential: 1e-150 (RND_MATEXPTOL)
+ rescaling of state vector 1: YES (RND_NORMALIZE1)
+ rescaling of state vector 2: YES (RND_NORMALIZE2)
+ precision of constants for MPF version: 80 decimal digits (RND_MPF_DIGITS)
+ precision of constants for Mathematica version: 40 decimal digits (RND_MATHPREC)
+*/
+
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <time.h>
+
+double nrand0 (void);
+void frandseed (long iseed);
+
+double
+lpsd_test1_random (void)
+{
+ static int first = 1;
+ static long double y[3];
+ static long double ynew[3];
+ double rvec[3];
+ long double eps[3];
+ int i;
+
+ if (first)
+ {
+ for (i = 0; i < 3; i++)
+ rvec[i] = nrand0 ();
+ y[0] =
+ 1.15332293276389647860084e+01L * rvec[0];
+ /* rms(y[0]) = 11.5332 */
+ y[1] =
+ 7.86361936577041429391864e+00L * rvec[1];
+ /* rms(y[1]) = 7.86362 */
+ y[2] =
+ -9.59832876597323340217198e-01L * rvec[0]+
+ 2.50191526820447899170527e+00L * rvec[2];
+ /* rms(y[2]) = 2.67971 */
+ first = 0;
+ }
+ else
+ {
+ for (i = 0; i < 3; i++)
+ rvec[i] = nrand0();
+ eps[0] =
+ 6.22795062914080487962406e-03L * rvec[0];
+ /* rms(eps[0]) = 0.00622795 */
+ eps[1] =
+ 3.18740533022684825004124e-01L * rvec[0]+
+ 9.37353183265263010305374e-02L * rvec[1];
+ /* rms(eps[1]) = 0.332238 */
+ eps[2] =
+ 1.34114066549657520897939e+00L * rvec[0]+
+ 1.37112680207817391021478e+00L * rvec[1]+
+ 9.09284711464199234257082e-01L * rvec[2];
+ /* rms(eps[2]) = 2.1226 */
+ ynew[0] = eps[0] + /* rms = 0.00622795 */
+ 9.99911280260497140622276e-01L * y[0] + /* rms = 11.5322 */
+ 4.56466267783293030389719e-02L * y[1] + /* rms = 0.358948 */
+ 4.76303657179166891018841e-03L * y[2] ; /* rms = 0.0127636 */
+ ynew[1] = eps[1] + /* rms = 0.332238 */
+ -5.54112616258510878069458e-03L * y[0] + /* rms = 0.0639071 */
+ 9.96821598309171222851810e-01L * y[1] + /* rms = 7.83863 */
+ 1.88399116837464845059769e-01L * y[2] ; /* rms = 0.504855 */
+ ynew[2] = eps[2] + /* rms = 2.1226 */
+ -3.92368890988536894184009e-02L * y[0] + /* rms = 0.452528 */
+ -2.28701078041458961515463e-02L * y[1] + /* rms = 0.179842 */
+ 5.25356147277256740562557e-01L * y[2] ; /* rms = 1.4078 */
+ for (i = 0; i < 3; i++)
+ y[i] = ynew[i];
+ }
+ return
+ 1.87990682693500129732031e+00L * y[0] + /* rms = 21.6814 */
+ 1.30213550872064624446121e+00L * y[1] + /* rms = 10.2395 */
+ 2.46871261992974920563735e+00L * y[2] ; /* rms = 6.61544 */
+}
+
+
+/* code generated by
+LISO analog circuit simulation 2.06b (Aug 28 2017) gerhard.heinzel@aei.mpg.de
+Input file lpsd_test2.fil, Output file lpsd_test2.out, Gnuplot file lpsd_test2.gnu
+ghh@hws393 Mon Mar 5 12:46:35 2018
+*/
+
+/*
+ zero 0.1
+ pole 100 3
+ pole 3000 0.7
+ zero 3000 10
+ pole 2000 0.7
+ pole 2000 0.7
+ pole 2000 0.7
+ factor 1e-3
+ rndgenalias 20k
+ freq log 1 10k 1000
+
+ LISO internal parameters:
+ internal MPF precision: 2048 bits (RND_MPF_PREC) [approx. 616 decimal digits]
+ convergence criterion for Matrix exponential: 1e-150 (RND_MATEXPTOL)
+ rescaling of state vector 1: YES (RND_NORMALIZE1)
+ rescaling of state vector 2: YES (RND_NORMALIZE2)
+ precision of constants for MPF version: 80 decimal digits (RND_MPF_DIGITS)
+ precision of constants for Mathematica version: 40 decimal digits (RND_MATHPREC)
+*/
+
+double
+lpsd_test2_random (void)
+{
+ static int first = 1;
+ static long double y[10];
+ static long double ynew[10];
+ double rvec[10];
+ long double eps[10];
+ int i;
+
+ if (first)
+ {
+ for (i = 0; i < 10; i++)
+ rvec[i] = nrand0 ();
+ y[0] =
+ 1.03417425760009383442729e+02L * rvec[0];
+ /* rms(y[0]) = 103.417 */
+ y[1] =
+ 5.49960697840505092774601e+01L * rvec[1];
+ /* rms(y[1]) = 54.9961 */
+ y[2] =
+ -2.22392902447038872852522e+01L * rvec[0]+
+ 3.81597301531777461688470e+01L * rvec[2];
+ /* rms(y[2]) = 44.1673 */
+ y[3] =
+ -7.95457894483044142182374e+00L * rvec[1]+
+ 3.26637284705751803829831e+01L * rvec[3];
+ /* rms(y[3]) = 33.6184 */
+ y[4] =
+ 3.16843598846870104359428e-01L * rvec[0]+
+ -1.28242087118526645285021e+01L * rvec[2]+
+ 1.48983635117700401584827e+01L * rvec[4];
+ /* rms(y[4]) = 19.6602 */
+ y[5] =
+ 7.13925936363737771639565e-01L * rvec[1]+
+ -9.32803024873575997255676e+00L * rvec[3]+
+ 6.92474395552374988608050e+00L * rvec[5];
+ /* rms(y[5]) = 11.6393 */
+ y[6] =
+ -1.41079946280815836171647e-02L * rvec[0]+
+ 1.76411538559131666157296e+00L * rvec[2]+
+ -4.09941918952344691377173e+00L * rvec[4]+
+ 2.49228067244039143945039e+00L * rvec[6];
+ /* rms(y[6]) = 5.11165 */
+ y[7] =
+ -4.63643951078101475109076e-02L * rvec[1]+
+ 9.69550099999588338527971e-01L * rvec[3]+
+ -1.60030456715483118846256e+00L * rvec[5]+
+ 9.11413893773142293781850e-01L * rvec[7];
+ /* rms(y[7]) = 2.08178 */
+ y[8] =
+ 6.22945857462875260502856e-04L * rvec[0]+
+ -1.24027844774927361142347e-01L * rvec[2]+
+ 4.61636778923473141831871e-01L * rvec[4]+
+ -6.85334702073203001693770e-01L * rvec[6]+
+ 4.33087180245228755075609e-01L * rvec[8];
+ /* rms(y[8]) = 0.941138 */
+ y[9] =
+ 4.06775537373803314227291e-03L * rvec[1]+
+ -1.23057208811098405605192e-01L * rvec[3]+
+ 3.65739331925435970013831e-01L * rvec[5]+
+ -6.05203599008239752941034e-01L * rvec[7]+
+ 6.29884019217132577216338e-01L * rvec[9];
+ /* rms(y[9]) = 0.954961 */
+ first = 0;
+ }
+ else
+ {
+ for (i = 0; i < 10; i++)
+ rvec[i] = nrand0();
+ eps[0] =
+ 1.09963864557172883519976e-08L * rvec[0];
+ /* rms(eps[0]) = 1.09964e-08 */
+ eps[1] =
+ 1.70997473272834082135666e-06L * rvec[0]+
+ 1.03692143869893427517672e-07L * rvec[1];
+ /* rms(eps[1]) = 1.71312e-06 */
+ eps[2] =
+ 1.78253746051597475328584e-04L * rvec[0]+
+ 2.32380677595372373222745e-05L * rvec[1]+
+ 1.58431751737622627485074e-06L * rvec[2];
+ /* rms(eps[2]) = 0.000179769 */
+ eps[3] =
+ 3.59813133547218049909647e-03L * rvec[0]+
+ 7.62196725361734776939539e-04L * rvec[1]+
+ 1.12826858599616888375909e-04L * rvec[2]+
+ 8.78842184705696600613203e-06L * rvec[3];
+ /* rms(eps[3]) = 0.00367971 */
+ eps[4] =
+ 2.67381910242847269433828e-02L * rvec[0]+
+ 8.26869471594200197800633e-03L * rvec[1]+
+ 2.01360904199803284576721e-03L * rvec[2]+
+ 3.45041110067685700144330e-04L * rvec[3]+
+ 3.13596941196716436649753e-05L * rvec[4];
+ /* rms(eps[4]) = 0.028062 */
+ eps[5] =
+ 1.27628163976407133110369e-01L * rvec[0]+
+ 5.49135386785089871164625e-02L * rvec[1]+
+ 1.98484934859710988748158e-02L * rvec[2]+
+ 5.69107101826107618059594e-03L * rvec[3]+
+ 1.15875047292676301571114e-03L * rvec[4]+
+ 1.26486945521247417203997e-04L * rvec[5];
+ /* rms(eps[5]) = 0.140471 */
+ eps[6] =
+ 2.76456036753196748373805e-01L * rvec[0]+
+ 1.63844860708100498076623e-01L * rvec[1]+
+ 8.45362863773692073233532e-02L * rvec[2]+
+ 3.68665394123939186542009e-02L * rvec[3]+
+ 1.28763989956705524912659e-02L * rvec[4]+
+ 3.23464258571102512348477e-03L * rvec[5]+
+ 4.42370286561799818810930e-04L * rvec[6];
+ /* rms(eps[6]) = 0.334597 */
+ eps[7] =
+ 2.73879355323971521520546e-01L * rvec[0]+
+ 2.37211874532050196111158e-01L * rvec[1]+
+ 1.78148428274873222526801e-01L * rvec[2]+
+ 1.15957373753373373794369e-01L * rvec[3]+
+ 6.41215950914302985406432e-02L * rvec[4]+
+ 2.87573235026066837288611e-02L * rvec[5]+
+ 9.45624128319538608299774e-03L * rvec[6]+
+ 1.73429939042556139874124e-03L * rvec[7];
+ /* rms(eps[7]) = 0.42602 */
+ eps[8] =
+ 1.23085854069711002632124e-02L * rvec[0]+
+ 9.34686518400146372025767e-02L * rvec[1]+
+ 1.44399590634281840340852e-01L * rvec[2]+
+ 1.61058887921429777516429e-01L * rvec[3]+
+ 1.47456168734661833779234e-01L * rvec[4]+
+ 1.13506617067785521288078e-01L * rvec[5]+
+ 7.18572879578562267969245e-02L * rvec[6]+
+ 3.44693368054946968380357e-02L * rvec[7]+
+ 9.65362229619917494276291e-03L * rvec[8];
+ /* rms(eps[8]) = 0.311048 */
+ eps[9] =
+ -2.91670247306267970006907e-01L * rvec[0]+
+ -2.55159448528591715766563e-01L * rvec[1]+
+ -1.72341547371169649232334e-01L * rvec[2]+
+ -5.79485593886401059725991e-02L * rvec[3]+
+ 6.79613344924523985049576e-02L * rvec[4]+
+ 1.83433626575543291836342e-01L * rvec[5]+
+ 2.67911575117128529308125e-01L * rvec[6]+
+ 3.04629928682971122934997e-01L * rvec[7]+
+ 2.80672610543033786601645e-01L * rvec[8]+
+ 1.75735805142254702358156e-01L * rvec[9];
+ /* rms(eps[9]) = 0.704088 */
+ ynew[0] = eps[0] + /* rms = 1.09964e-08 */
+ 9.99999999989765830389858e-01L * y[0] + /* rms = 103.417 */
+ 5.85448303147780718959120e-02L * y[1] + /* rms = 3.21974 */
+ 2.25369821241975815303325e-03L * y[2] + /* rms = 0.0995398 */
+ 2.57907804956185268284736e-04L * y[3] + /* rms = 0.00867044 */
+ 5.03964020282466854368257e-05L * y[4] + /* rms = 0.000990801 */
+ 9.81023918210109789840666e-06L * y[5] + /* rms = 0.000114185 */
+ 2.57309749087191638205109e-06L * y[6] + /* rms = 1.31528e-05 */
+ 7.43644874155278312735896e-07L * y[7] + /* rms = 1.54811e-06 */
+ 2.02392982154030339665238e-07L * y[8] + /* rms = 1.9048e-07 */
+ 2.61754473467469363553839e-08L * y[9] ; /* rms = 2.49965e-08 */
+ ynew[1] = eps[1] + /* rms = 1.71312e-06 */
+ -1.68198677749980277453802e-09L * y[0] + /* rms = 1.73947e-07 */
+ 9.99999998124022169874095e-01L * y[1] + /* rms = 54.9961 */
+ 7.69905045333711544963394e-02L * y[2] + /* rms = 3.40046 */
+ 1.32158973558155890704524e-02L * y[3] + /* rms = 0.444297 */
+ 3.44322501473490584970214e-03L * y[4] + /* rms = 0.0676944 */
+ 8.37748344272230346471635e-04L * y[5] + /* rms = 0.00975083 */
+ 2.63505164313549931114686e-04L * y[6] + /* rms = 0.00134695 */
+ 8.85489272238472508158737e-05L * y[7] + /* rms = 0.00018434 */
+ 2.72003077561341858164931e-05L * y[8] + /* rms = 2.55992e-05 */
+ 3.83683073814904877682951e-06L * y[9] ; /* rms = 3.66402e-06 */
+ ynew[2] = eps[2] + /* rms = 0.000179769 */
+ -1.87478919894376624622559e-07L * y[0] + /* rms = 1.93886e-05 */
+ -2.09239982599206658867025e-07L * y[1] + /* rms = 1.15074e-05 */
+ 9.99999037954570059091600e-01L * y[2] + /* rms = 44.1673 */
+ 3.43310347891692251012971e-01L * y[3] + /* rms = 11.5415 */
+ 1.34162661624854490221722e-01L * y[4] + /* rms = 2.63766 */
+ 4.35133554149551351274250e-02L * y[5] + /* rms = 0.506466 */
+ 1.70878981129488790619294e-02L * y[6] + /* rms = 0.0873474 */
+ 6.85514275160114284729647e-03L * y[7] + /* rms = 0.0142709 */
+ 2.41617234854039262864161e-03L * y[8] + /* rms = 0.00227395 */
+ 3.75324681274365778089154e-04L * y[9] ; /* rms = 0.00035842 */
+ ynew[3] = eps[3] + /* rms = 0.00367971 */
+ -4.11276534124143521449322e-06L * y[0] + /* rms = 0.000425332 */
+ -4.59405513299382065397536e-06L * y[1] + /* rms = 0.000252655 */
+ -2.11103671796321695962402e-05L * y[2] + /* rms = 0.000932388 */
+ 9.99942378819654214123270e-01L * y[3] + /* rms = 33.6164 */
+ 7.81439853087134473342356e-01L * y[4] + /* rms = 15.3632 */
+ 3.79952687137442521777665e-01L * y[5] + /* rms = 4.42239 */
+ 1.98473432503483102128042e-01L * y[6] + /* rms = 1.01453 */
+ 9.87016838890094808431582e-02L * y[7] + /* rms = 0.205476 */
+ 4.08087905854408957608840e-02L * y[8] + /* rms = 0.0384067 */
+ 7.06967405056280260015198e-03L * y[9] ; /* rms = 0.00675126 */
+ ynew[4] = eps[4] + /* rms = 0.028062 */
+ -3.40265302718047872593669e-05L * y[0] + /* rms = 0.00351894 */
+ -3.80519770089554275646831e-05L * y[1] + /* rms = 0.00209271 */
+ -1.74718695314966633764990e-04L * y[2] + /* rms = 0.00771685 */
+ -4.78035672132966189883839e-04L * y[3] + /* rms = 0.0160708 */
+ 9.98455494313708245731161e-01L * y[4] + /* rms = 19.6298 */
+ 9.69495252964178393084254e-01L * y[5] + /* rms = 11.2843 */
+ 7.55823862808095504826297e-01L * y[6] + /* rms = 3.86351 */
+ 4.94210553110644742456929e-01L * y[7] + /* rms = 1.02884 */
+ 2.47436485910456842529965e-01L * y[8] + /* rms = 0.232872 */
+ 4.86617623601759649863203e-02L * y[9] ; /* rms = 0.0464701 */
+ ynew[5] = eps[5] + /* rms = 0.140471 */
+ -1.88064922021257489386039e-04L * y[0] + /* rms = 0.0194492 */
+ -2.10657483300338225226917e-04L * y[1] + /* rms = 0.0115853 */
+ -9.66180309680497668973801e-04L * y[2] + /* rms = 0.0426736 */
+ -2.65248070625893357559299e-03L * y[3] + /* rms = 0.0891721 */
+ -8.60187370064416057493878e-03L * y[4] + /* rms = 0.169114 */
+ 9.78081304490830575663911e-01L * y[5] + /* rms = 11.3842 */
+ 1.51371631662000148409259e+00L * y[6] + /* rms = 7.73759 */
+ 1.45117965063678327927021e+00L * y[7] + /* rms = 3.02104 */
+ 9.26791298750758446808194e-01L * y[8] + /* rms = 0.872238 */
+ 2.12805397554288695654208e-01L * y[9] ; /* rms = 0.203221 */
+ ynew[6] = eps[6] + /* rms = 0.334597 */
+ -5.08142240947043395281091e-04L * y[0] + /* rms = 0.0525508 */
+ -5.70645196872719451945184e-04L * y[1] + /* rms = 0.0313832 */
+ -2.61273445899179266581295e-03L * y[2] + /* rms = 0.115397 */
+ -7.21091474413940775774300e-03L * y[3] + /* rms = 0.242419 */
+ -2.35206394945239000658234e-02L * y[4] + /* rms = 0.46242 */
+ -6.03647298611043096610507e-02L * y[5] + /* rms = 0.702605 */
+ 8.28316537290167089099942e-01L * y[6] + /* rms = 4.23407 */
+ 1.59879449090784422722329e+00L * y[7] + /* rms = 3.32835 */
+ 1.45977742224345071123660e+00L * y[8] + /* rms = 1.37385 */
+ 4.13975528263623848236135e-01L * y[9] ; /* rms = 0.39533 */
+ ynew[7] = eps[7] + /* rms = 0.42602 */
+ -7.62639006307564650251477e-04L * y[0] + /* rms = 0.0788702 */
+ -8.60525281366654206035637e-04L * y[1] + /* rms = 0.0473255 */
+ -3.92735453879055573911272e-03L * y[2] + /* rms = 0.173461 */
+ -1.09456774515701846376425e-02L * y[3] + /* rms = 0.367976 */
+ -3.60856932828152869538449e-02L * y[4] + /* rms = 0.709451 */
+ -9.38329650683611704459635e-02L * y[5] + /* rms = 1.09215 */
+ -2.71344607384770184186407e-01L * y[6] + /* rms = 1.38702 */
+ 2.83178245818690403394832e-01L * y[7] + /* rms = 0.589516 */
+ 8.01743167780795391991907e-01L * y[8] + /* rms = 0.754551 */
+ 3.46815631487601873174976e-01L * y[9] ; /* rms = 0.331195 */
+ ynew[8] = eps[8] + /* rms = 0.311048 */
+ -5.63060291907531205530437e-04L * y[0] + /* rms = 0.0582302 */
+ -6.43843021371736296074378e-04L * y[1] + /* rms = 0.0354088 */
+ -2.91229346986437781789258e-03L * y[2] + /* rms = 0.128628 */
+ -8.33882331184381660556425e-03L * y[3] + /* rms = 0.280338 */
+ -2.82966566981860351091115e-02L * y[4] + /* rms = 0.556317 */
+ -7.61586921402223412613324e-02L * y[5] + /* rms = 0.886436 */
+ -2.29749051885100855333936e-01L * y[6] + /* rms = 1.1744 */
+ -6.41607860395749695547441e-01L * y[7] + /* rms = 1.33569 */
+ -4.64310221765327115144999e-01L * y[8] + /* rms = 0.43698 */
+ -4.08952897236250812109629e-02L * y[9] ; /* rms = 0.0390534 */
+ ynew[9] = eps[9] + /* rms = 0.704088 */
+ 9.41771813006329633714985e-05L * y[0] + /* rms = 0.00973956 */
+ 8.42870661187725399378889e-05L * y[1] + /* rms = 0.00463546 */
+ 4.51982523805608457937101e-04L * y[2] + /* rms = 0.0199629 */
+ 6.85382010343668709301712e-04L * y[3] + /* rms = 0.0230414 */
+ 1.25464229662007344354339e-04L * y[4] + /* rms = 0.00246665 */
+ -6.65839096109422030787845e-03L * y[5] + /* rms = 0.0774992 */
+ -4.56927215863728720774663e-02L * y[6] + /* rms = 0.233565 */
+ -2.19879910305822107530734e-01L * y[7] + /* rms = 0.457743 */
+ -7.85068586700891144510576e-01L * y[8] + /* rms = 0.738858 */
+ -3.62931337954710155224358e-01L * y[9] ; /* rms = 0.346585 */
+ for (i = 0; i < 10; i++)
+ y[i] = ynew[i];
+ }
+ return
+ 2.09861132274782616561919e-04L * y[0] + /* rms = 0.0217033 */
+ 3.91085881911858824284933e-01L * y[1] + /* rms = 21.5082 */
+ 3.19581417187890788195496e-03L * y[2] + /* rms = 0.141151 */
+ 1.16373972646361867320508e-02L * y[3] ; /* rms = 0.39123 */
+}
+
+int
+main (int argc, char **argv)
+{
+ int i, n;
+ if (argc != 2 || sscanf (argv[1], "%d", &n) != 1)
+ n = 20;
+ frandseed(time(0));
+printf("# fsamp = 20000 Hz\n");
+ for (i = 0; i < n; i++)
+ {
+ printf ("%23.16g %23.16g\n", lpsd_test1_random (), lpsd_test2_random ());
+ if ((i % 1000) == 0)
+ fprintf(stderr,"%d ",i);
+ }
+ return 0;
+}
+
diff --git a/example/lpsd_test.pdf b/example/lpsd_test.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..cdc5c1a9efb1dffa2554f4596b64ba67fb41b182
Binary files /dev/null and b/example/lpsd_test.pdf differ
diff --git a/example/random.c b/example/random.c
new file mode 100644
index 0000000000000000000000000000000000000000..9e740c43bbad02a4ce7b67217868f92836c1ecf1
--- /dev/null
+++ b/example/random.c
@@ -0,0 +1,320 @@
+#include <time.h>
+#include <math.h>
+
+#include "random.h"
+
+static int iset_nrand0 = 0;
+static int iset_nrand = 0;
+
+double
+nrand0 (void)
+/* normal distribution random number */
+{
+ static double gset;
+ double fac, rsq, v1, v2;
+
+ if (iset_nrand0 == 0)
+ {
+ do
+ {
+ v1 = 2.0 * frand () - 1.0;
+ v2 = 2.0 * frand () - 1.0;
+ rsq = v1 * v1 + v2 * v2;
+ }
+ while (rsq >= 1.0 || rsq < 1e-8);
+ fac = sqrt (-2.0 * log (rsq) / rsq);
+ gset = v1 * fac;
+ iset_nrand0 = 1;
+ return (v2 * fac);
+ }
+ else
+ {
+ iset_nrand0 = 0;
+ return gset;
+ }
+}
+
+double
+nrand (double mean, double sigma)
+/* normal distribution random number */
+{
+ static double gset;
+ double fac, rsq, v1, v2;
+
+ if (iset_nrand == 0)
+ {
+ do
+ {
+ v1 = 2.0 * frand () - 1.0;
+ v2 = 2.0 * frand () - 1.0;
+ rsq = v1 * v1 + v2 * v2;
+ }
+ while (rsq >= 1.0 || rsq < 1e-8);
+ fac = sqrt (-2.0 * log (rsq) / rsq);
+ gset = v1 * fac;
+ iset_nrand = 1;
+ return (v2 * fac) * sigma + mean;
+ }
+ else
+ {
+ iset_nrand = 0;
+ return gset * sigma + mean;
+ }
+}
+
+int
+irand (int max)
+/* integer random number 0...(max-1) */
+{
+ int result;
+
+ while ((result = (int) (frand () * (double) max)) >= max);
+ return result;
+}
+
+/* A C-program for MT19937: Real number version([0,1)-interval) */
+/* (1999/10/28) */
+/* genrand() generates one pseudorandom real number (double) */
+/* which is uniformly distributed on [0,1)-interval, for each */
+/* call. sgenrand(seed) sets initial values to the working area */
+/* of 624 words. Before genrand(), sgenrand(seed) must be */
+/* called once. (seed is any 32-bit integer.) */
+/* Integer generator is obtained by modifying two lines. */
+/* Coded by Takuji Nishimura, considering the suggestions by */
+/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */
+
+/* This library is free software; you can redistribute it and/or */
+/* modify it under the terms of the GNU Library General Public */
+/* License as published by the Free Software Foundation; either */
+/* version 2 of the License, or (at your option) any later */
+/* version. */
+/* This library is distributed in the hope that it will be useful, */
+/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
+/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */
+/* See the GNU Library General Public License for more details. */
+/* You should have received a copy of the GNU Library General */
+/* Public License along with this library; if not, write to the */
+/* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA */
+/* 02111-1307 USA */
+
+/* Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. */
+/* When you use this, send an email to: matumoto@math.keio.ac.jp */
+/* with an appropriate reference to your work. */
+
+/* REFERENCE */
+/* M. Matsumoto and T. Nishimura, */
+/* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */
+/* Pseudo-Random Number Generator", */
+/* ACM Transactions on Modeling and Computer Simulation, */
+/* Vol. 8, No. 1, January 1998, pp 3--30. */
+
+/* Period parameters */
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0df /* constant vector a */
+#define UPPER_MASK 0x80000000 /* most significant w-r bits */
+#define LOWER_MASK 0x7fffffff /* least significant r bits */
+
+/* Tempering parameters */
+#define TEMPERING_MASK_B 0x9d2c5680
+#define TEMPERING_MASK_C 0xefc60000
+#define TEMPERING_SHIFT_U(y) (y >> 11)
+#define TEMPERING_SHIFT_S(y) (y << 7)
+#define TEMPERING_SHIFT_T(y) (y << 15)
+#define TEMPERING_SHIFT_L(y) (y >> 18)
+
+static unsigned long mt[N]; /* the array for the state vector */
+static int mti = N + 1; /* mti==N+1 means mt[N] is not initialized */
+
+/* Initializing the array with a seed */
+void
+frandseed (long iseed)
+{
+ int i;
+ unsigned long seed;
+
+ iset_nrand0 = 0;
+ iset_nrand = 0;
+ if (iseed == 0)
+ iseed = (unsigned) ((long) time (NULL));
+ seed = (unsigned) iseed;
+
+
+ for (i = 0; i < N; i++)
+ {
+ mt[i] = seed & 0xffff0000;
+ seed = 69069 * seed + 1;
+ mt[i] |= (seed & 0xffff0000) >> 16;
+ seed = 69069 * seed + 1;
+ }
+ mti = N;
+}
+
+/* Initialization by "sgenrand()" is an example. Theoretically, */
+/* there are 2^19937-1 possible states as an intial state. */
+/* This function allows to choose any of 2^19937-1 ones. */
+/* Essential bits in "seed_array[]" is following 19937 bits: */
+/* (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1]. */
+/* (seed_array[0]&LOWER_MASK) is discarded. */
+/* Theoretically, */
+/* (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1] */
+/* can take any values except all zeros. */
+
+void
+lsgenrand (long unsigned int *seed_array)
+ /* the length of seed_array[] must be at least N */
+{
+ int i;
+
+ for (i = 0; i < N; i++)
+ mt[i] = seed_array[i];
+ mti = N;
+}
+
+double /* generating reals */
+ /* unsigned long *//* for integer generation */
+frand (void)
+{
+ unsigned long y;
+ static unsigned long mag01[2] = { 0x0, MATRIX_A };
+ volatile int kk;
+ double res;
+ /* mag01[x] = x * MATRIX_A for x=0,1 */
+
+ if (mti >= N)
+ { /* generate N words at one time */
+
+ if (mti == N + 1) /* if sgenrand() has not been called, */
+ frandseed (4357); /* a default initial seed is used */
+
+ for (kk = 0; kk < N - M; kk++)
+ {
+/* printf("kk=%d kk+M=%d\n",kk,kk+M); */
+ y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
+ }
+ for (; kk < N - 1; kk++)
+ {
+ y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
+ }
+ y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
+ mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1];
+
+ mti = 0;
+ }
+
+ y = mt[mti++];
+ y ^= TEMPERING_SHIFT_U (y);
+ y ^= TEMPERING_SHIFT_S (y) & TEMPERING_MASK_B;
+ y ^= TEMPERING_SHIFT_T (y) & TEMPERING_MASK_C;
+ y ^= TEMPERING_SHIFT_L (y);
+ res = (double) y *2.3283064365386963e-10;
+/* printf("=== rnd=%.20f\n",res); */
+ return (res); /* reals: [0,1)-interval */
+ /* return y; *//* for integer generation */
+}
+
+#define DFLEN 64
+static double df[DFLEN] = {
+ -1.8370939043417808943e-06,
+ -6.1196945589982073985e-06,
+ -7.5587701579763848592e-07,
+ 2.5562408350808373654e-05,
+ 4.0758326159745981627e-05,
+ -2.5131872908960921123e-05,
+ -0.00014876214901936348016,
+ -0.00011893811742986104748,
+ 0.00021896415862725504128,
+ 0.00052444196762459233553,
+ 0.00011727581338660759027,
+ -0.00095051536688676254964,
+ -0.0012527267103039769756,
+ 0.00048357364652680834428,
+ 0.0027983445801144007873,
+ 0.0019952329898457937768,
+ -0.0028299913322717383303,
+ -0.0061757318063923029715,
+ -0.001405108411686069795,
+ 0.0086145480140459759784,
+ 0.010554659468493140902,
+ -0.0033774389783623680766,
+ -0.019616438026572809159,
+ -0.013502662547377131163,
+ 0.017384385505533751742,
+ 0.037703687212232158688,
+ 0.0091398768675829813617,
+ -0.05266831278756242539,
+ -0.07033376546222222311,
+ 0.024396819308122194912,
+ 0.2055912251957120361,
+ 0.35282470939550079266,
+ 0.35282470939550079266,
+ 0.2055912251957120361,
+ 0.024396819308122194912,
+ -0.07033376546222222311,
+ -0.05266831278756242539,
+ 0.0091398768675829813617,
+ 0.037703687212232158688,
+ 0.017384385505533751742,
+ -0.013502662547377131163,
+ -0.019616438026572809159,
+ -0.0033774389783623680766,
+ 0.010554659468493140902,
+ 0.0086145480140459759784,
+ -0.001405108411686069795,
+ -0.0061757318063923029715,
+ -0.0028299913322717383303,
+ 0.0019952329898457937768,
+ 0.0027983445801144007873,
+ 0.00048357364652680834428,
+ -0.0012527267103039769756,
+ -0.00095051536688676254964,
+ 0.00011727581338660759027,
+ 0.00052444196762459233553,
+ 0.00021896415862725504128,
+ -0.00011893811742986104748,
+ -0.00014876214901936348016,
+ -2.5131872908960921123e-05,
+ 4.0758326159745981627e-05,
+ 2.5562408350808373654e-05,
+ -7.5587701579763848592e-07,
+ -6.1196945589982073985e-06,
+ -1.8370939043417808943e-06
+};
+
+double
+filrand (void)
+/*
+bandlimited white noise
+PSD (fsamp=1Hz) :
+1 / sqrt(Hz) for 0<f<0.15
+dropping for 0.15<f<0.25
+<1e-6/sqrt(Hz) for f>0.25
+GHH AEI 23.10.2002
+*/
+{
+ static int fill = 0, pos;
+ int i;
+ double ret;
+ static double save[DFLEN];
+
+ if (fill == 0)
+ { /* first call */
+ for (i = 0; i < DFLEN; i++)
+ save[i] = nrand0 ();
+ pos = 0;
+ fill = 1;
+ }
+ else
+ {
+ save[pos] = nrand0 ();
+ pos = (pos + 1) % DFLEN;
+ }
+ ret = 0.;
+ for (i = 0; i < DFLEN; i++)
+ {
+ ret += save[(i + pos) % DFLEN] * df[i];
+ }
+ return ret * M_SQRT1_2;
+}
diff --git a/example/random.h b/example/random.h
new file mode 100644
index 0000000000000000000000000000000000000000..fcff032db882a784a43362aa51d40b586a8aa23d
--- /dev/null
+++ b/example/random.h
@@ -0,0 +1,7 @@
+double nrand0 (void);
+double nrand (double mean, double sigma);
+int irand (int max);
+void frandseed (long iseed);
+void lsgenrand (long unsigned int *seed_array);
+double frand (void);
+double filrand(void);