Skip to content
Snippets Groups Projects
Commit af306525 authored by Axel Schnitger's avatar Axel Schnitger
Browse files

Add an example and update documentation.

parent 682d1728
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
.gitignore
+ 4
− 0
View file @ af306525
......@@ -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
Makefile
+ 10
− 0
View file @ af306525
......@@ -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
README.md
+ 82
− 42
View file @ af306525
# 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
frequencies equally spaced on a logarithmic axis.
......@@ -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)
* [Usage](#usage)
* [Example](#example)
* [Options](#options)
* [Documentation](#documentation)
* [Theoretical Background](#theoretical-background)
* [File Overview](#file-overview)
* [Options](#options)
* [FAQ](#faq)
<!-- vim-markdown-toc -->
## Install
......@@ -25,12 +34,77 @@ 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-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
### Theoretical Background
......@@ -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
```
example/Makefile 0 → 100644
+ 9
− 0
View file @ af306525
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
example/lpsd_test.c 0 → 100644
+ 407
− 0
View file @ af306525
/*
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;
}
example/lpsd_test.pdf 0 → 100644
+ 0
− 0
View file @ af306525
File added
example/random.c 0 → 100644
+ 320
− 0
View file @ af306525
#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;
}
example/random.h 0 → 100644
+ 7
− 0
View file @ af306525
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);
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint