From af306525cee3df256b8ee73f129d5c3942fad868 Mon Sep 17 00:00:00 2001
From: Axel Schnitger <axel.schnitger@aei.mpg.de>
Date: Wed, 7 Mar 2018 09:49:10 +0100
Subject: [PATCH] Add an example and update documentation.
---
.gitignore | 4 +
Makefile | 10 ++
README.md | 124 ++++++++-----
example/Makefile | 9 +
example/lpsd_test.c | 407 ++++++++++++++++++++++++++++++++++++++++++
example/lpsd_test.pdf | Bin 0 -> 23142 bytes
example/random.c | 320 +++++++++++++++++++++++++++++++++
example/random.h | 7 +
8 files changed, 839 insertions(+), 42 deletions(-)
create mode 100644 example/Makefile
create mode 100644 example/lpsd_test.c
create mode 100644 example/lpsd_test.pdf
create mode 100644 example/random.c
create mode 100644 example/random.h
diff --git a/.gitignore b/.gitignore
index 204659a..14bc376 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 87be461..ab1e7af 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 a367fcb..8cc6a60 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 0000000..1095f4d
--- /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 0000000..0c84a92
--- /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
GIT binary patch
literal 23142
zcmY!laB<T$)HCH$-THRjZ!Tj61BLvgEG`=x1%02?y!4U`1rr4Wg&-~k1qFS#%$$<c
zA_aZ7oWzn;m(=9^lvFM|JFeoAqSVA(u8KKJYrW6w?kv5#r|ziaHnzVjPdqA?QrVcO
z!hCzyNhy|1IR|zWDOk1%{{JPiQ*Mpb-Cd<3jsm|P?>nlrb7%Jc9pdNZxBWH!_r1LS
z>#wK(Z^i%rlYjf))cF5}LI0oq`M*DB|G%@@AmJUd@o&G{$JbZ1{`}GZ*xqo;_5Aw%
zQ~!y6*1qyjG^qOTztG3wx0^Zd$JPJ;72Wsgczxa9>G8R_|L(^9lvy=p%k=AoWhY;H
zJpH=*e*XSH%jLhze=lAT|M&7Y=Nt7OSJ&S>6@Kgc;^Z9jqW@m1Tc^ai@B8xqE5H3e
zr8fCf`gNaI@4xx)_WEl~`^|qA>{EzS=Vwy;Keumb-PFIgCbm59d#U}r;4R0Y;^RL=
ztM-?%Jly&CU+J6ef{qn$X50?n`BsPXyzJ{A=a=5ku`k!YZuhjA@8#;Rw^<y2xO_Qh
z+4@&sVDIG4)%LFX|1B3yzo)8EtKs<ahgF`D_o3qB<&GZmT!#YpSMBPO<2)~W`g)yk
znb17zW$Pbr+pJ<!ZnJJ~AHP-0`l_j+TjuY6n}45;bACr0_w%~Htc$K+T59jA|JQPn
z_<dE4It|B%A#crRsy^KLc(+0eKdYil{L^J$Ju2QwJo`MoPr;`A&iyYsoc*%j%ciSr
znZNt0d9%Zz;^Uvb^-lQE$<byz|I{=62?sas`n2wUpS?53@7Wq(&#SClzWZR>?~S)-
z&wW3C=Fe;U^{lO;`(y8I+g!VCw(k8<o!k9)@9pXQHq-cG-)_|ekV5WD>*P)}A2?c9
zwsK|L|M*7feT(898oTS;4qP=)Ul>>0`o84$`?Oham*>6z_ptq4;z{%KoM(^w<?oKP
zyM0#Z-Kwc?kLTsb9cjP!Wanga?gQVFx0_#Aj<@ZfFS_%x-IEHQ&XddgtdqDusw<Ue
zeNil)&*@iZ{I&A`32wiqc_}yR-h4}&ZC~>4QtQ^2GuIx!ohLK@?Xj2Nm%h&1@b|{s
zv%l(Ju+Fo)f3{^ah?i=yPVVse&NCAp%t(4+T$cMb-Syddff+y2W+eYSQQrGe&SUf6
z7wmg}r^xS0dHQbcl>RwW&0kB~ZQi&1d|CO*zVkM#=I<z&Us>>bxAWJDzHiDez7s$B
zFVc4M{#z$E?pVJ5!)GtHvbUZoT508;i#r13m$dO$-FqW0R}tO6KKdHl+2aR~CCg=f
ztCU&c<G6=M)v0bnYt4nj8}EOd`_`g-t=!|&wMYL<?JvlC@w;x1-LDqwFK=gU<@DEi
z)I4kJFWZfpZ;sb|lqw9I*LG#H)MZ<FcHZ|%+co-|+jc+L`TX;O_d**r&NtsnJ!C7-
z-dk`(uX)xd3yt|-7QFgtlQAoHU31qb8;yA}&6}QoUQqht!|NF;HvE?IG0KbXKW<Qd
zeldTw$<sm)yO)MnKUW4+Y`^a4T6syks)+a0y!)4RCQBzCQ@9_fEC2GB!J~z*&)lp_
zxmI&-v)g)ix2f|v>V8IiI6eEpbbr?RBbMb84W3qBs1oP=S|jsii~IJ;9Kv$#N1Ej(
z{+-9ZX8Yp4`}@V7R-UN(sIO-$&z}3R*50>LW{Q7Rz}Ff%ld~1CtyXB9Z+>*rRHNFC
zKWce&rEJKQ=Z+trNKLVoXMg*9_bIcdf$Kkb@~@hA|EjibUGJ3>AO7rL9Bvn+Q&|81
zjo8!9<f@>_<>%G=j($kX_o<LEy7Q*@mo1;f#7#SUGEP07P;)7H)$_*-j_;1qIrVtK
z@mRe(pDQnXKc5#kk9}`_*}I=M8P%7wx<A=u?9Md(UL#|2)O>zHpd9<t$#NmK^V!cD
zGZrEoG@pH~aej$x$P{(f#p>4=y!mO%D%}5hVSsb?{cwZ%nKLHcEBUm@W>b%U)wfqa
zPi*=%`TT;qKay9rpI=b-OLEJ6kCfjLO1EdbZ1Y#$bk+RvK8^Fu*Jc-%pFi$sZR9iE
z5fLm~?y%{|^M>g$PxWVeneTXLkq!T~e8;A0JKm}{7JQM-dZ#Wgs0wP{1xf<<pDwb=
z6t&&SmU{ku*3{y~pEE7x4=v=sddl>x*_H2G`qT?86zc!n^;a`g`|x;%r<MIPqfhyN
zEw<)sE$^$If8_s<-~WI8D*f?v|KErG|DVg({doQVPka2IG@ldquG^^0ls$a%RoRny
z$F%Eyys!UrSN`AM^8deX|Ns5#YyAJ`Z$3Tl*dJ{hXa6pK?bS=?S^q!&b~yk4_wDun
z|L$M6<kHEhe_s^t`M-MZt!G;}pD#N6+Ww!*hRMs=`m1NV_be3s_jiU^ab2}aTw8?Y
zL0gW;OPdXhk1KA=oRilwW4q^!O@e)Irq43?6sFD?dAT5=I<Lj>o<h#zxr^M+=9rz~
zlDtw7b)w<9c;n||;f7xE=jw006X|Qb+&^L0vz;95n|@o=NhNOF^(Ild)p+Le3`O_C
z!@2vI<ReeCt8L?We(Bkf3ElFlZj$dKSmqfiytV6_`)TjV%(=<QQLc+0h@=Rdin*N=
z;@8z{8TDl4b&Z3w_RMMNowCz*;_;cy2H|@-0(dvy(3@m%e7DfK{o6bljVo7`D?hn7
zXF1E_tA-l49WP5|-*tGnMPv5k3c<|_G{kNklTP%Co1!tV#ipWcv58OjricSaq!rI~
zrfuVTR%*8W=CeJTUC+}!?&%z2d>1J6bKyOXBr(^GDQi;%PuQM%=J~+R>C?Hgvl-`n
zbvBDlath~~zTsAis_Q$Ig3#`QcVFiSEnOF~w%e#o=)@gYiQV-qonZ^tafq#*Y&2a^
z#WPi7h4WV7#gP}EX-V(c%`x|W=?0OK0^^Fy)=RFo>}sCrSba|T`(%c=xo_Jd3r}ZV
zII`1OUGZFHq4hblSu;ZCZl1aMP+rJ}UpbTKi_W&?nB#n9%gMb`C--WGiN^ol*gah|
z>gA22JGtIka?H8h82aW>sPFA=&kHwadEYoIs?ykhH^y^9qioZ;EH;)I3BheFZ=P{&
zeV1o<z|n#;#rfE*vds$}(q(u54A`M(_UR>m#o5mshGpwFoMSt_s`=>c1sSat)_r%i
zoqU9?!$Y2TC4D#BZrMG#pyxdQ6s9h<o)yMr`NIFY&Q3lOnEhGh_k>j^g4SM?x!mnL
zIZNZ873Y&AX8Yx`snzFp^;+gE>6z|s!zb7k+H^r-s>jB~wk*nF{TBP@<pz8>P%~rF
zon=g&p{o`;=Fer(+u?Zs*rSgbnn@QMHpHaQmfC#q#m^~K0&3CS(<%}ikBY5P-C$Zh
z#pJ0{+FL#DiAmoyGQHftd9mHzY_{&8Ylr9W5}^meUJZVSo)jC+H>jBNOKP)LM5Ix9
zK#`!Zba80={Rjn}Z#7XznHKMzV8%Ut@!iJ2Go@xf)90pt6Is%BXi7%*H7V73%T_qH
zEb4V}HZ7gvn;~#2L-yW-Pe*qAxVKT^^_##+W;b$p?|wL{>vUhg|9h>}9BIYGYg6X)
z@jkiL?D$~U3R9~sPa>`zGrN9t?fVONF8*Whe!D1kS)Zz#b&KtxbZbq4OMKOGlb@T~
z_O8%+aH0R=gTpuG&p5iH&Vl>b8^#ChN7lTx;^cKYzxJ)?eAcB4)$3MVp8Imb?g(2(
zk7fHG&s&<g#p4HqdSm<n&xTck(MDgFN_0=!zf{p$tvoB<aG^)D_sRKFm{k22zH@yR
za(RW(@|kCu!{Rru__lG;yanB=Zd<zl-<0;w`BBik{<V8%WdG^UmsLGpImBO?|FL7;
zZ<P$E5<lbG#XBx5Ubx45s8LI4O`*WDLft#>99L|9lkv{?4s(;-%8QOOu5J00EZdUl
z8yaiIF}XkFbK!CJ!flqCzgG43$hTPL&YTBQQXQ~sl~F~r!P0*(Vty}p5`Js%x}TOR
zksmI~UH0C+oOitftK@t~$^X-uem88?+`(@*S;R)8q^L{I%vmMqi+8|=sf&J#TwBGe
z$9mt@D@oj9_9T_k05<EOMJ4?S<>!=t&i!B({LTC0%Oz*@RnB$2$`vrG?76l5O>Fb#
z3DZ&@Bp1ECtKd^6qVsVH%LZ%D2)6)RwWu$b80Kq-p4C2ZVEL-bRtDa=37g!j*Dc9f
zn^;r3(yQU{5s4*gft-aVO!Ip;NPFrXUvMVWb287NCwZ5rY}EH$rOKAWlk}{23e%~F
z?rCcizI!js4*tI6d72vc^P6@v=a-+Cik?}1Ugnzh`=qjl(DS?Zon{a>@$bu9(zEnQ
ztoOn2dFB0GOWgV7#Pl+zwMn<^oNifRpK!DI!;VsuEf+Gc6o0UK*pcR@dwa>5sjS@$
zcecyea<2IKCSNUKd(jsescSR-3bNk#E~BY(K_s%OOyJp}MG=NF4}I8^XTP=LoD*`#
zeWSLz$VTo+voe8at5&RZuQvKM;k%7Sap0%2q?YxwCLWhOyRkp=R;>Wnzqxa|OG1s$
zmMxK=$nR-2XLe1+x#W+QAJ1=`vt+&2rNFZ%y&GcB2WfpdT_n`6d@0JZl2yV&f2;7v
ze`|LAaJ2Z)5!PCN%9Cr7UHZ0VQrnbI3))5(d!F=enVk6Irox%?48?6)r#|o<GrM&D
z+|AQ1rzdXSxA@_n=bjICJ-nO#`d~z&#il(k9Z%M@OrBA9a;dHRlji~pbJ7-ZojtiO
zXRp-5S<iCR9JfBnQ7c{R?b9@uIZ?Uw^`e3Urn0S*BiLkT=N#_y_;T%BanzGH8;b*{
zs2#g4tR=y}WWk*oK4A+BUFs60MNUaymJny!b1CQX!J2cu20LRGEOeZiowP!v|M9Qy
zF}@eqo&ED-t>sjokL;f7HE&1A{t(%uwAn-CPH%ugdfwZEPp-61+Pdb<vqs6)Cshky
z_wS6&;-9`}&!#z#DnyKaI(f`6W__#?bxrBaoe({tsf)8ByB*w|gZi}AZQ}Kw$$f=a
zBc?TFW#i6CbB=|sF^kj=e8nN1Kcy+sb63JoF~x&|XQ%d_v@l)3A=-XiB#+lgQX`Xd
z+jF-Wmzk5<uD|T<G~z6NE&7Z9(3JbG78bUr;ulp{swSKgIyfbd*Kg7x_WQYl+dS^;
zYO_hX-=FhDNoIy{tG;9N9gYe9ySzB#Uo5R;c<NZ#5?z)(IXIYMd!>lnuSGsg>VBM0
z5`XqyIN;6o;KNpn)vk)`6`EFZ)*ULlzDKjcvp}>nW7f0VcJhB$NjGXQWOx>PCx+$3
zip4@s{&_8$EEV&<MfaZcw4Pbeu+2L4!!}Q@=iOn+>)ao{$?#v-wSo6k%tYJkJmx0a
zX^V`Da^)qOmDsN*?GBcbXPj=d=aO{e>qTdHFaJArR?J>~U6;`GbY|fN_Pb`q2(@-U
z(>33{e2Z=5&jTLCi7GWqww?2x<b1b6*}?Ke?DC`&eX<{0@7>BMIvLs{SA5gq<;f|I
z8BsG?n*XhkDvr(a5i3ufn&dW5v;4iS?ZtOz-KP{N3miGRSO2$J)QuOLgg+jQ)&1hS
zzGbJ<@nfy~GQy4Y#i9yQn%y@i^ko$Z-r41J{_?@nL@&+qw;?n4Ek8V)>Hp^9;8TTv
zrSF7%s}gvisJYWtBVlFP0+}er@&yMz>sU*tnyKY7x4y94wB>-S4&QQyf4jG4#}(|p
zIFGNma8GE$+-_ku);^~~hOmTbr~S+hIZV5Iv0&w#tv@;b#?L9%xNtPpTzPev)4Q#f
z>KoS;iYCd}oq5F2-|oih?|jH>!|oV?jbiQ9tkE$|Kbw3tZO(Rz=A{-|M3i?pU2F_c
z30r329rXN3H}f&8|AOleI~@I(%q!5gNO+Z_<_y!wpvBW1^+I+`a&=E{J`ifc^~WhO
zs6ne~_oD^ES5LbBd}xtqaa1wGHm&T)^(LM}PjW<(CrIe`1!~_p`AA@^5{K5sTU(UQ
zY4SvA&p9b{mM7xRqAlMQQqDx}IT29rT9~W5Ug+$ML#Cbjtl2{R50s}h^K$H2e2?q<
z3Jd<DPP2<<uorN59eTj^J32;1pmI%!3|n}_;-UzN1O=|gkF35{Se|M6ck9%{ug<Ds
z>?cJnS1yv6>F#sXu&9-<k!fad*6wMpW+pX{R5uz<e{y2y!(?r~Intj4UDXwviyjE(
z{k15Y<$GqKZ}<X>33(F^6`G&@!`1oaTdZlbf#h3*U5-p$i=VOOFEH?yzulM;QmAs5
zZ=Khw1&zy9EAO%=Uf4V_?&!{g)6%@QTXe~>`)xjX(Xju~jG0HDA5|=NzCJ0pNNUm3
zM;0^F{f-tEI%h=5h@0<VK3Z6)eZ{n^N?K{<qZQkI`?pz~$o(2BbMeuEl_LLmZ67Oa
zJ1D5Z=OJBlG-yxGj&xVUjOP0xT}SsGoT|yESotPI#=vaHwpPKU=Bs*r0=I39booBq
zGMRn!VoWoqPS-svoruCU$DGwgbsK^Vaz8CRY%@hNw)Dr=qaQ=smGt@gtZHvRN?75X
zE!J*sD>Uu=TX}V##V1|MU#t$&{2;5f!s5e=g0$Yp@-3RWeOgu}4uw1y`Hxt+Iy(DD
zaeUwRxLT*|r?uUNB=-I1&o4jl_}vGNH;c<XvNle>H?dE%?4<Rd>q_D>>K?ZhI`|91
zyZ^rQn}2|>g#AAE-e-K%rlrkKlu_*6CJ}6yy8Zc&B%x(tZ)D4RHj7psmpw05eRyvs
z?;|OhtI`h-OwkR9T`eU4@bt`rH(!LVNc3F2c{bzu^|EaLeW#q|S~TT~I{f<*V#5!K
z^skG#xZq6d$#l7vlbO5DZhlt0f4lnfr)K<2J}(!^Y@TSSxnqGila5p1rFhnPDmnKV
zm(STWo89DTiC6TU{`vQA9sG2B$10osUN2JZ?|+#3sczpt|IkV2Cy8FQ?Nca`<Esz*
zeWPch`;DtlWI|UgS9`zG>9&bwjHuHd)42B`%NONZiLQ2TnH{sI%YVm7*)>08X8$>U
zDyiRX`sPQEGd~4tyX<&*W|33S>{HEJPgY!HzsHuClQH>vjNo;XpEpWZ-B4S;$G0%I
z=(F^f%-c(!++43&^1#72)n-n{?T_ChLP8asQzG7arA%lIZ4{1-6U{XB$PCrE)|Rw(
z)ve6SuUD_T>1cUsb<i<s8Pj;l8xKlCYb$xJw!Ye&TQX_YzQ1f!qgL8V+04?aJ!~{(
zmB9XkKUak0-QJO&ew_KnEwj*^rTcP}o#U4JN^i-Y`Y+dDw?@b<2MOz-sUle$CWeYF
zble(tb3$9fYNel77ettgcYl7hKuh}ChkIJe*PpGnoNoCz>f`e#3(saZZ(3_~fR}6Q
zB?ohZ{iha~XS8cI-?zG?7%J!Fd?0k{@^=ACr~0|Ow?sW%`sL%^a7)A5`GH<jPhUQ~
z`jq&qfJIaH1t0$&{xqhvW^L%Zi1ZI3Pd98_x8&6E*B!f7%{kDwVD%~QH+wFrYiqsl
z&@uQ`xYIPl?Bb!>5AKCm?v`Cw<$Fq|bjQ6?s|$-wBenkTQa-n;=y+v!=+yo9EEj9Z
z&p)bS`tZG`SbG}V^-ZfbtehU6y5ZvUX(3{U=f6E!ku;re`C6?jFAEN=>RMxMv*y%`
zmEjfqq04(@C&h?PD6KFH{TS`OZll(o*FX3oeRujE3gG8bKOPg(c{@>k^{b^;AzrB_
zw<lg-8Tf6hb93l)|F<2BR!Iem#C+f6JTvV2s&{7|hF@Rx?!x`6Q&(*6dlDMywYmR)
znE2dxmp6r^%5HO<aVBJulECG_*w9_d9jh99cPdQp+0>A-TJPOW!RKM&o$qdF-1ju)
zn0#>8rxbyT?5SRD-%O&`INlb1zr(4N@3rSDZHID>ic-!?^CX&gzbQUnWE2s#HNDbq
zMb9lkT~XCXY>jaowkHn!n)E2@{N5X1A9_o=m7Wy|oT6?2eV6Y4+?a*#`oWXee429f
zj+TGqp0%HhD(fusc5oNC9KEAvZoIkj^Bnd@J2mZHWzQns_VaViP49I0zo^>&-<ry!
zOpDJ7v;-e<ZMofLz`{@%V8qce+hqBJ8=b{m+GUTdoCD4{#(kP6dpuh4i0G`jye%tk
zy%Size$A�-}qHwmXMdcueGJPYTH^OHg?_`NDxL(bkQN&q{;_r>KP5xVi=1yzIln
z8T(0Gh&kWSWc>!)4M9gI{!>s{!fkLeH(*NOi&vK?=m}(Z8*$xv%fS6|%~1sfrIM}D
z2?bBBt_MtAc~LGrC2Py!sWYV*q?|$*xoyrmbtNOrXU8rrGn3RT6YXQZ#(S4&+n#=?
zX|${8Va|(F?U`M%-_MEOxL6o?YF|*=;?C+<J02JAU^yBvZ-=(3*w+gG6RQgMCkMFo
z{FZGFnttcK#FSl6%8oBMlj^!qi#vPnDxNi&rAEa9Gj?CGS^hmLW|K63%G?;)(?Vg-
zGuRg@M6and)G*FaQhVBFx_-%>4bP@dk-oNPlF^-6;bmR3_0MV_?e&$KdR^PoZk9(u
z@%hLlK_Awf=iTviqrumd%fht<2W-<@ck*)`Gu5=QJdyM1j>FTNI)&dEk14$1XnCJw
z$Mn#;vf{vgf7{2*(=X~>cb;`Xseb;49&Po?p2eDF^EGpGJ`1$V^B?54TfyaG_hXK>
zfbGN0n+jgMxWDY_4w<*C2|hIv5*<+=8<P}F+}bDHUh?70#NEDg`p-KsZ}hN=;9Ya|
z>&+a7O|OftNpqNeKDlw-zK@k94wYS}X7)Efe^+wl(8Oz+e|GfFOn*1)QS`<G3Llot
z3iesUv~g-+Y}Z<z#`bmFC0N9^XLMCwZ9T*L#I;G!OtNLis_-=H^-(i6bOoFE&T_o|
zgXh%RB8?3mH$vT;4}8+`UA^;5r@(~8D`mVtXRNBd<nwFQJ*ft-E4CBE?|zM#CDP5l
zN&V;}#>c7=CT})vIk4(srPRKO0Rf!uaZ6*57R^a%+TZ=kQiw0)5NGj9DW<8#%NHow
z?JM5bKiMbHv~yNUh_8~`PEM;;OiqdJn?J7VG-Y9%|H<ymo&!s?lp;<9L}*R5Ynu`>
zbw#9Fm}=Idm#15&O5ZA-8h7s7lWiCJShl1^3o%7)*7aPZlQMhKl+eE`BR5StU|=t%
zna*yY5mc8t<I0moVXtO}#BR3wxas<$YZGFwc-c<YI;}0cLFL+pTPYhJS=`Q0$yk}Q
z+EGAnRk3-{vW^=88x~mH4p5o0)=KH*Z!HeKnR;up+gcx;<#UYSdb8{Cfm8RNaHt&S
zQEQ5DOx$X0)>Grf=NNJF^{aHp##nh<j>T8*Ok<cmKZip@?_`^T^3<c62bYGrep6HV
zsD0DuiL>kN6cO*l*k`go?y5>R`OG)FDYQYy^rX;}$S}QgzO9dpiWw$Pe#|{JOO0LP
zO5wU_{*dl}sZ5U3BVXt}C|YIqDdd`G)qK;ulUegw(@$p24?FtxTjH$+Z>M%GH(mYN
zRQ1R9lNPeKuBH5t(k*oQutD;FPKZOoUC&K@%exa47#S}%Z+LKC<D>bZ-a4O4kqQEr
zPD;sr6*%cNZ^BCx^Hz?RJB6|<f0cW7WGQL1>^<5(Uq*9*K*H<^T>BVKsYWioyw&&G
z3!4i6KE4-$%byeo#8qq(n{QXMiEsJ!`AZJ7uJr3;E8yF^cT%EbkKd$HGsCO<ZqJzU
zic{mEcSn$AOM<{MtDBG0H%wczOz@U!@}<*1qJw8#+OV|sb6fBPu{3F=Y3nw;KKY^C
zr*s#?l85h?w00buK1uA}%njndKe<Gz|Ni75p0hs1)<&4;WR!$@^2rOyX{M9*md?{y
zs#$Q`{NOrqKaUl!UGC4B@j7Lx_nA^7p$d16d14>R#rJqR+=>-<=YA;7lYeM|OXBI9
zY#g^aWtZLLQ!UX_T6XwX$bvVlcY7}w9gC5ide=t9F#OB~j%m52nh|B2w{TCtd*i{{
zYfkLTI|UEh2~5gWkC^ag>7SO9(>@7n+}`CkBWH1O>&dWhEEC@33T)baMm2J;;p94r
zof<ydw%%N`&~npSMWMg#rr~QBD^B>DvT3VG;LY}_LV-7LCx%ZxkaPc-Zo;=iUp|@~
z-G5VQ!MBZDcqZOD?qm63$Mb@=g};>Y_ibQKTXFrM?7eHYvZ6(u-o_OM+upS(xK-ah
z@6fSro%fBN%W*qA4s3{Tw?ABV^gDC=-Uq*qB<8o<O3Rh#R_rOb&gb8{yX;c%j%;~J
z>B_e^f(~SR?`gQsJ>T!aR^uLCL5D3<^MxJst_f7=P85B2>srK)?6!#8OmWs7Lhp7}
zcsyE{Yumba*%i~UBW~%Z81~k=Mq4Yb7r#}lv|hGDpq4i;=f2pVs)f5Bs=KiB9g6?G
zPvYVF+F)+2#BJIuZgb?km1Ml(y@lhzx3>yR3GVluj0J9(X$d#hwuk@M+_SYx;7+A;
z0=H+Mq4=YyfX50NI~N=Y<evDL@5I;CuT3$t?OS4IH!>IHeq$=$8!LFLUEuc}<-~R3
zSH24*#RtSEt=r5Rw;-FVE8x1A%EEQJSxvFgwq<FD78Nr+-t$K3LQAaogVtTv4c
zwg}$XCt0@FPIJ+^>+bGXzW&nb+!goTkD+|;$H!`^pKi~$Z@eDpcsd}_`rnUJ`(7Lt
zcpckcdN*{@_maEQ7oB~->!qv5?bz>ECCuKuu4|cmx6-lq!P!QsnXRr)0lPhlt}oSE
z*CNt#tF`Y?fOiVdBC{-uQ$7-=KF#lpoEulfy6m{TB;9P;!pkb@QHh&Ng&TfFFg#NE
z?Y^X9(of;l6~FxrCrxX3<*ZtEBsl)S^F0oyC05#U&s%4?S>;>7UTOXVe9BD>8>Xon
z&Z_)->sZF6iwTEn1FYsu)MEJdWUKb+jth1%Hf|l?oOXoAWk2yUP};)(SQ#X^$F=3#
zLob~_&#oGJYd_*E@4K|%zO3N1+1tBH*RMRSD02J3(!BOlpE|GD#kA<J_;${%kzd<G
zFY(qX^ZX;3{9OeHS+aKDl6A<PzD;<=+ecp+ugv?_rt$XW(#`22BKBO9YuO)hFWFmg
zq9<3(#9oF^z0CGdK@#I`>81s&ht!PXC%7EiCN5B^vR&po$JV(q2Ys)U6&);FU{#i&
zdSjcj{=#kI|K@MVJ)Yl`yPVZrE$_CYLg8&~ACa8f6Q4BXTAyRdsg9nwqTu$)(hDuQ
z>fyG++umk`Z+LrpE6b5w>lC%T-$nk6iMRW9v)+jl7iZXETg-fGYF|i$%0_L`2J<6-
zr-~Vshf6tFytUlEWx=)&agz(TaQvC1-O+j5ua)<i*^66T+vdio$jsz5J~w5<+Z0i@
z8)Zc&_4TR?h1zc9re}2AoyVfrb<$w^8=>93HW%Ey9P>`T$h*7hK)J!+r(MxXTd!~B
za>&Y`K3nvWPXT96-T8XXm3!AkJt%eUWdAmM<L%e)SRZW{?v8(%u5x{|(4M!eC6xRt
zCWk6rynS-%ij7)VHJ(@-C$I%G?wPx*dt37EJ9iwueA~5?@yFY#cSUNVU-~cFzWsD8
zUsL}3oHIv0zvC3#nfqemg4-9iJOAh0og3=)X4|DbD?WYkPvzW|8zJDd@@>(@<`1_U
zd6zCSRK3SPeQ8|ly=R;+wmrCg>;uQX+!=ijZnwSU*qd8hEpyvz1NVX4*luB?$GI#I
z#P+y-c)Q|J)Z^`(^8~AIo7-{h)lK0!aQpRx#v5<n7QJkKaQoUGrh@3tLI-rs>mF_A
zpK2|2Aooib!=19*C2nOiG+on;9&y_z-Fd9=VcSLZFJ3jd%@^*lP3PD%x91^iRkVj-
z)AriKO$ToOIKFN-_mPHmzeOtICp6so_KW4ovNqA4+X{To4%P%ZD#_%<x-it3zZdg(
zSd@9d;JLWXz6p1a_-TahxbWdnDYuMb#CB8l;A!F?ejKU{F4kMHeZg(J)1s$tdF1O}
zJib$4(v)r5B9hM&d+YBP>~pJ?UhI8%daQG{%nt8=t<2jOYsl=_x~T8)hxL0R^B#J)
zGH>5KgX_TTs077|%|9yV{gtkZE&ZYYrB3uiKy-NMC5MWdP}gR*Oc$M-woeXrFYL2f
zt+|^$lj&pTgmqrLUg?+e^|L1QS@G&FSzVX%@D$@CY2JBVx!hl-mplL1B$p`L<8b!(
zs^*z*7Z}aZJ5nO>z*nDF>%&Fw!?K^)=5GIx+0c|P&Jbz#%lklPSnbW;N0qr7FCXbn
z>#|RoCoEpK-6q0pQ<q0U!%{V_I~$hIVYR>2^NVd-hx&)5zdReuuB7~(lNQ-xY8bU^
z1y{fs!=^OHLy_znC1s(Nizj^gkjj1WjPU!nOkNN7#<ptiDBqg$NubB_$J^N-@2%4n
zeYU&Jdw%5X?2}2n-m#lp+81An6O4WPDt6(%P(GQ{cG=yA3i}pEF7!F@eH+)j)*oez
zTLtc~JKE>+<JYw!#RIF?&lTag)B2*RF+AH}EGhKx9*(k7&wnnB{(lZUEIWIQpGTs?
zwCC;vmcxfSH-FfF-ud6Q_3Ucrocy@F-^_~iJCOSM`l<JeE%nu7DvD<u;a_5z%@kAd
zxz_cE!mQq@1rKvKa=Pue-0kosK6qQh)mdlr^iNyg-sRA9;6qwk!_ivD(+`h+Yd!sp
zU-ssFsRJMK78n?0)~*%)w#2$h*=>z*h@E}oyjoeY?E2ll+JU}zC)g?-{T(H0w=r*Z
zP~|G=vbQViTnpnjyxO^Y^S4<o2B#OW8XT)~4>Z0iX1_D?aF48_UCqBu=C_+wH#g0A
z|75+o!E?iF#}9pvJovi(Sd~+OLH*Cbl7}lN$_G2Xt`V%ia;oRidGlL`m~{<SS+GX7
zCCv|=zU0FbwIi%wj%}BH_I8z1AA9^uttS;7!S=$*hNsvTPUXBbqqS9rd0x&3(K#)B
z{#_Q4f1G4q^f8(){FoKyw8{3KOwCFi#wQvZp3Lx={L%cI*`ytM`#bh87kv9j$W?C6
zU$L&3I|>#>rtHm&jYYTKzsI(_KvuA?ukKg=@3b2g1y->;K0dH6YL(_c^6%Y~4etcA
z_mt@y?A3c<(e1onQmroW$C1wV9S8O?darC~WjDI=YnsEy$Tyr~!kX(Ioy%|QYTsew
z(R7qOv+?l30BPStN1O9|#NupTv|Tyd@$o_P4n}GIeeKdaK1L?5<m5lD_~TR&)A8+*
z&Ys&owkOw2-Mip%;?l#Mvi#b=&lYdimEoSo|8(iC6ZviZ?DE@%?E7xH2Om?~zTurU
zYfDX-@uYyzQypU0c5Gu<mv1qvL}cQVhZF89PqJCVqPy@(;?ihNweZY`-_B-ynjqHv
zFzrmx!JBQ&GGV7CHGF#Ts3+IyRP5{Tz7sxeFlsz9rG4rt<yAESElkTU3&u1n$5c9G
z^(^@G;FhA|Y3E%^CpR8DwaYW<bfE0?tJ$Y+c)b+Jwd=lrL(cj^&ndPYl?9=5K0Hf|
zjq*HF8S#_d>0I+Cc858Or}=OcoQRGnH)1Mz_h@6j+gx_>ZlSV$3zc|gJ3Xi@P}G_D
zJTZ3THRDy=Eq3;(9*Pj`l)JXqqh%3$JeP|`RYdr+k{;0#>D~h`6Jx(Ql~idsz9_Iz
zOS~KDR8e(6R?%g7d(s){2QL>N(%{-3ck!CZtnAKxk_+UnsdCRP_ijDS%NCKg<BmK-
zdi+z{`#nMY;w!8JHyqrc_V9}H-J7$TQf;OOyjru=FvvFKzBKciss-oumQ)$I_y>Gd
z@L90n)xmR?8(%#x=lOCVW^Qbu;#GIm&JC}BdFp(0$teAI?cz71xLUF52S<;uxN;!=
z&AN_y#*6>ot&dA~wg11_QvOfeh5LE``5*rJzKtt2FQpi=Of5Ap1+>5ow5kn6D<~*f
zT3DEX6o9x!@bz!t6>qK;CGNo`i6yBDHZ~vu&?>k@0|n6PIA=!%Lk0chL<K_yJ3EjZ
zgw%J-&nrPvou+G`YoK7NU;tS|XJ-dmF$Z;n8Qd}j{gBFnR0aKj#Pn2<S^5EqMX7lu
z3WhL!`o5_tnTbyM6$;S`1_}lW#%2ZzMuuhzv0&BC`FSO&c_qaP#&GrSMfs%#NbZ1`
z4;BS!0eK`?K|iD@F|W7)WLI(}OxRPQq^K0^3Xme_U<G}b)UwRv)F5{!kP@gjgHnt0
zON)|IixogY7KG$0b0l9umk?UOorP@)A;=ntLPG^ZQv+i&kYUi+s+jY(s<U1sblZo$
z+gI&c_5D??ul~yt^(BX=wF-1BVdCsgyzgT$aRVo(feMSj1XYdhmQWFq7zbBY*33x_
zz6Tf!(_JJIXB^}W2`N!KbjDd`dddFc)zU4M``&;5*RWpe`MT$M&*yEm;!ti0Ss`e6
z!aVa?%dUR=x+-?3d#gDn7EXS4_G9kK-vx3`E1Ly_N_KD0dv)gbrSks`0mqbDcAU5(
z$^EW2^6z$LgTrzX!A7~KT288St;=^=;Vri1>kKyDV>4&X4VcgKbxWr4^E2@)ubylA
z<s49-acqsJut4Av#VesLzi!PGes+E1`z*P}Uy}pAUp+f7===A*U*jGAa^!4c7Tgud
z7k5Lq|JOF58qb67C&Ny*SbP&|4za1`RBm7JKu_=)$IHGc725Bb-XwmkN_2l1@G<J)
zgg&Pbmeh~o9=pmvOaEz8P?6a9B3!=sO?apCeg*kz9>pF#2Hv`deT=12D)P#iygM)7
zT$;jsQ)@$Y#2Oc_J_9elY!17HP8WX|o_xXoq|EsSd-7Z<ze}%*EpM$^d{*XLp^&>p
z`z~JbN1_{+FMPCg)h~k_-b+XJb>*<$2)M9UdYw~fV*P;y0Us)DtDZcc5FgcjzNhX@
z_qyU6zklm|+Gq6ONW~2)DdGQZrU(AHOljiVqE*9pQMNZeSIzFL{A2Nnq6at@^v6$8
z+td~MMDeThukR&uCWYVjT$#Dysh_+3^^D-9OD`E-GMzhr-Mo}pi8J4=DAQH=uK8@G
z_^WVDx#!zE-ruxWdM@((iu~)je*C9r+D5Hi>!^8AS<^76%WKh{7p5i)q}jhJa&Fsp
z@$)$;udk~LcGv7H3%L?^#m_&ie^r40uLOq48E+;|&fO}ds}#r*blF+$;;hMc4=c1W
zt=yAz=8LiQGt=`gqpP3I+{`U!x25{xqPx4N>+IUxHRqttzU1>irbYZ&9HH~=*uCTN
zNBxh=?Jg4aKfLw+^)*Y%!+RDlnRbT7?&Kty(j|S0&L0Gv3K#^{!`V6bw&le>;ybcF
z@Vep#^N{}4K9zOehA&$dvugw&c5Gwl)pd9*GNp6Y(l0OhVh(vv6LRm&pHsLc-^zIA
ze4Q&&>c5|H6yMcunfk8yP~_^y{$u*vmi=b3vlmfuZgaoHcU@@hwKs>-j~%`hXyhAI
zBqu)om0@<F$oDBaFB~>FPFbmyq|&7%Xr8w6R)PPMgRGNQoiMv+c6!^DW49KmtW}+}
zB)59fw(CMy-mch{Vjxzsc2UYSPM(;Sp8>mkyjkxvYvgX3A?lExBJf%B?#7#5C+AM*
z*5G<JU99@rf|tqN$9+_ex-Z*ScBj1R=ktx;{%iN``F^?YPvYzLbN7C}`O6)j<o+^i
zH}A`9&wu-?g&wc_vteZ|r{3Z25AUv*%l`j*dDV2SMgI2t^5vJxPWyTJ{N4te7uK3f
z&+JrKs+rpT()F_Rtm0Yyu}zm!uBTc_rkbyty^Ht96vm)VwW&uljOJc=dDd!j<nrro
zC;d(?39i(cl-J((qNM6%ldjlyLodGltI|J8b_NOGUdr&gFq`-B(a1+Zh0Q0rdTW|v
zmo1Z3-kqWTt=I0c;GFg^6~>ny%{#k`d4`ZU<K<QL$x-<t`M2tu^Vc7oy(C}sgqmu8
zp!L%Aj1D>C`fG)Dwd_3lXg0@+y5=WeLzn%D&N#Jv_8&2O<}Gew#{5(Gdn?awDm})$
zw)*<Lt!ByJtBUqLzM6Pi-mb0H^7XpCw;QF8E;G%y+r73jV&c?4#r`>Qsc(MGzhAlI
zXHt{+T#k)O*V%q~bg%QY?Rh#Yb5`!GyR*KV+(}5=@GWqfVUm(UZ_`7|sk0tVbX$CN
zUP6TN;SD`vg1yt8&OH3YK*(X@y9ps3sf?_112zU0y<aV}d{WS!ovSp?P2aL-iqwgx
z(tj)exY;gT#Cv9TcZWw!iz0hT;GVE$m+e<IggyMQwP4|ots2<|0UNJ#^Rg~D-owhe
zdRI%N!i2c>?7zOO*0FEd^6~tvUe;R~f7UEHr<o9BEpCwG<&`Kgk1HYkz{E3Bp0Nuu
zW#Sf0bYh(qV*X?<_aP%TjaN2vl22@8;pt(VtSfu^WN`nUzh&7v{QdRwEPuS7d;a)5
zMV*^>3}c_|-}rHTr1ptCiMbPvEY^lbuZmrn>djj&`E2KmiPt2P)>dvRU02pvxpdpf
zJ-5E@2|M~zLTh8zR=3s3l`<Ey9D3j6EQ~3WxE-keVD6r<-eZSu-9F*>gu&2nbHyE*
z*Dsftn)Q4?`$<r~AkO329sxH+Ll-A$>)VcpX3M76^*`Q}_k4%+pJ^f)k@IdJZAoBu
z&E>S4%)0wunA_ve`A>JsZSrFce4g3W@cP)cDs%t%UGHxmwQN4?R_a`mntL{?c=i9?
z;m_vJsC)l}-Rzj&{A$}i`7QaeqLm33jh`()`RC|lxrlvt${);Kd3N$~`<q8U-};rj
zQGOqrkq8t2Cx>AEUv951&VIet;*4yl%(MG~_liF^n~2W3xa;Mtle=1vwJy7QO!?UQ
zIq^H!S(&<DT|cYCyQEulk$cD=!#dHA3LluBoca-Z@?~IWlT+cLw<f#q`W<|m($GCO
z&+WCK=Eru)4U8^-8m=Wjo50Qze<^uUpy`cs2V{f29v@DT{G=#&WyP_A1sy%=n@<EK
zG@DIo+f(5@Y0s1$2bP6|H$C;w@{Vp=wYntskX1s6j<?aS#JN>WGgqgv>;4O{i+Cuv
z{qppzpjBa!|E||0aocTwdurX<@Ae;5>{5PxpMO3^{>P(_>Dso{-^IlJthPNlyZ+tW
ze<vpUE(?CX|914hOB3hv`@YcooGz`h`hw?%K*y7QTob!Q`rJ2P2-?EG(Rx8W`|mT4
zKix~MdvQ-!UH|x{vQJMut;)FXIIF#2;(O{L<r8u2u;RhvNdiyBRb~Ap=D%qR(X0wT
z6&x`A`Ri3jZYHWZ_Z`f?Sv05S#MZxm;;R`q3YC1lztebybeE^UUy|BJ-ZL`Oc$Yml
zWPMk6Mxn02Yvp#8we`O=mNl*S*{Pkp)pPnn-lY>w)NF4!&zIKVoK`0LTv}Q}c+SFg
zSN*Qv&7PMOzaVZ!?234ubH3Re>r!T2S-sWwSLvg!E$q#+w?u!5)|r<u%Ww6wKkr2P
zy&s6`IEz@!n8uaMWzaX*++ySAyIoe@H`?1DXT`LeC$UDH?htwN=3|89xyz9+xh&s@
zv<0Z&IXUw!zhogtkCExSyXtp+dw*|Tr8xc8(QgYK8WhdXsEDdJWPf?Y`+>cR>-417
zRR<UuHhvPApjvhF<-#pr0yOl^75&ZZT-M)D+7=l(TPJS8`O8jI*JmByzq{g^_u{tj
zZ+&Ma0*y;g?w6~5w|r{M?7H{wP8hS=e*F3ASgUq@<*T_j>T{3Qi~is6Zrwg!)uqAA
z+c+QodA{sSRc+DP-=|Bmuj%iuj_^D8rRN^|v$M8W4W36w-nZx7_s(KRvTcR@>7^6S
zcmL*e&~EB<{PJRA0{71;5ih?i>o|G(+nu7D*H=%PqxXB~Uoj6`32}W57QY16u!kG0
z-?}Ue%M6JUUA!UUY}n+ffjcB-%S)F|xWT$zK01`mo~ymkVh-ndqZn47$+??v9f&m+
zmgQvmd((3xgZ$zDCjVYOE%erUX2KKu<HJeOjemDMc=&wxbiT@*V;A;r{Jzu7`_zNZ
zjluc*B)Z?uo3MAg`}a%v@jreY_u+r~)O4d%QhV{I70aKb&NNU{TF#+-?}Wf1?YP99
z3u^Y*9tyskC-XikabZNo>7PCa&M_MBWPjN8`FXPGpN_NkbMN~+cl(&{!T7<n*-e&7
z{PE_rV#XVr&#TPe`M~gk@tLwTdFL6G9J3jv4{>)W&N$r>!)AAzS5M*gCV^>7vYeMZ
z(b0-ITFM#TxN5<=9E+TTmJ;vpNobcB{QoiYTtnBR=`;W4$Y0@|KmA>L@Y$7nr<7}X
zUpltT&og#O%=(yl-<n=NUGumnEVKAl?%h+H&fT;Smvru!(Hr{BYhh62%2grZSz24S
zI<4Il7P@-f)hk)i&TC&*o$Ae;TeYM&ICO=2yIT4g&NH6bKGiN?m!t(BOOBkq`)$bE
zYc|vO8824-5_V#i)A@`qB0a*Fnx#Zve`(cOv&*~7*eh<*ojCW2$upAoB+q|5#d%?Q
zh3sY-ahWyOHy(79X+J5E!uYJ(Q^-~4pznjrGZb`~au4uk7)BRo>Q$Ds#q{hBS$4)l
zQ8SPG{I1<A_A0E;<U5&izSQO9-Lo}OY6<&Z=~>L}G@g62Bt*m3c)Qi_2s6!7d1usa
zDjl=n{dwujN^t=_r^5|TW$zY$7H*2RZ@E<SEkID`(br!F9r-(Ix9J@Cnjmuc;i_LZ
zLwXchIscutVOs3{a;`1oxjCXn9c~v@jJiVBo^PEZxre)WuHwfT3twGljnNS0*>l13
z=^@FNUQ?CM?3xhnQ<(Uma>u7bE=?<+Pn3+?mLtg(%xd^sOkCe~%DvWxx)cB9rI_OQ
zJy&!D1a7YT`nUVpb5{Av>em|{&6~;JANA*!g|o}$fXyqeWxahBRrh7%M$>$KyHdIP
z59ircd!3YJtTw2=uxTN)SW3lV=J)fn<RdMO>>v0UuDn##B_7&$yYBkD7JCbWu){Gd
zal#X&{{~jY_#RTZd~f!YiL6qt6PKO{e0Jl-&4&+*e0Q!>^6y-E#PpQODbtfCCqv(6
zewh>c=je~EkD||(mu(dGyxOA6s^V|l@#T7tvgwmO&ys&nzm_i>zq?Mg>P*&IzwNWG
zYo5sq>n$;h<>TsllxZvWb<)15`zzAF7V~7f{|x@aYIkmP+`bO^W5#=wj~XAJZ+%CO
zkALZPPO<I#0^aO!+k2$!hq1_sW81pgTE5<x(BYjDn|=DUp+m-}=#)uQf8LmJJt;hG
zZse!1znf+M*?yLeyJ4RHu6Ti2%RT2p@iXr`Y=0d3HUGC^wg0@D$8$G7zvp>GUc12~
zpShOjh2zTdaAC1XnS|n-b#FAj#g*}`4v=D5y5Pzz9xvS_mNjMC6>A){QjfVYa<8-#
zHk%mYEO2_^TGi$5??PAk*Ss*<)e&KG#o~c)wNb_U9V%h5vuio|-1<Wc4*zw1RK2>S
z`~H&lIlfcV`&{j0^{WNBx)YCXuFdLtw{N-s{e64i8&==AVEmQ2dhZOo1uLq`jz{kL
zzvah%54G*0Hv=9kE@67f-BL8Q#Q11HMoB`8X_RQ%gwT+-uUf5F&#w5Iq~&!qDmkF}
z@t%NHI%}W3{^3@qd|v6^@%Imu{}}!=S^wa+hwi4wa=+jE$~?TKx^b>lulJX-TPwYm
zZcjdOQt-aF#;yguMMWCt1sZmroNilSS*Eh6xZvRt$;i{wx64>0aHUtwyQ%kW<|mz*
zTf{yaC%wK_W@Dss+VlJ-4YQ_@W;wBovXiu8`3=`i6n8P~Dvr8PHkV;%%7g8%-#PZ2
z_#x4DGr1yb?Vjy@uOIACEz!Tr#V6>TSJ(eXW`V*5fshleXXSY$%TsPPeX(XzShG>k
z{lor^?|%1N)mg?~-1788wqbscRR7!r@8E;SR$MfZYL%IM$#+`Jhgnlr$*#JrwdfH?
zAa{aHt<CxsSyz&D1J!$0)%<MeUS(pZ`^`}ye(l=JoAzqIS-^K->6G5wt6!Zg|1|oF
zg_KU!72-SXC+5I#Y3}ip=a2W7YwXz8Cw^kb{NIzO^(h+lrzsqnlj`izE<b6Dz_oqq
z0?OyiZWXm={9*RwWW96z@$#*YFU<{j;bJr;_9#=#0;!hYlm8u3cxtdws9$YyYvHMf
zk9cf4eV7gx+RszIvp%LdoMZJ=`-j;_o(eeam@E(`#&Mpz{^%YL)6E$tRU?lvsq!AX
z<k&bV&P!A<WR-$i8qcq|)xshvOyX?%;iit=TN3AXxoq4QDp29aC>3|gsb*DKJD>90
zC9AKQzFPOn;;zh+sHTs0E$xl!((J{rE_4_F>AvTa)p$rzsr0}rZhglE9KW_Q-&h;)
zhs)xfHuKd-SN5!n`m=0BwXu(<1mE>=#Z?W81~wNs+gY+&g;E4HCN3!Z6Y3@Br7S6U
zJ^X=|<rR?wH(tMdIqBxp@K{!_gRgcud)iOBnQ?i)>eQ_bGZ^PwbLsG0d}fZ#8jc^I
z|A?LWe&S=a(TQ_ypKchW6&-)QS?tQA{rO2ON0<G!dbczF=aKpkTh*FVI$Tx?85(J6
zrmx;+JfmovTyFZii?8=sExMijr*@witHaMmTh67=8$NF6SduI$WB;}9OY7;r9_b=2
zGp6_>^W`nKEX};rzQ3gAOYNI@L&s<HzeMonO|4xXbMow^(?XI$ZZ}puu<(gXoT<?N
zO=G?Udv@G!+5Mcq9p^Du*4>fHwy2j|-+eVFY)O7wWqITK%}mci7a3^^-zjvwaUgk~
zJCFBYf0YmZPHR@H9a)(+)$@I@>t6q?T`!OO+<dxBN^sgjdEXT^7j_)Z|9E;s_s=ia
z*Ly5Uwz$5{{91PY@x2M>Pd%C*$^Rmo`^WD2o|E^U$W=_9p|#G-(#^JY=DD*$eu3Od
zELBRL=JTG{nZ=iXy4$-~=I-62#}v;cedkDCn!KdCsc>0wQ}w>2eNU}+=7~=feZAG}
zc0u;mt>;cf-2K)2h;L$t((9tt6)t}jRqy369Zs9PCFI)F2a8=C{>_bh^`6Cfr^)mn
z8-a_}8?Vn&UiL94D(#)G^7EY!H?*{u#7JpfYZt!9^*3S4pSX}4f$<u*nUg<V)>}8f
zH}~ZG0u6O($uoQ6ziDiU{};79;a%qGU;i#_->7}#&&l`=eE+K=ukJ~@bgQrb1c&eW
zjjP4q9!t(WzcooeZ3A2S5xefcD!;yM^*ua0{M+nrPKUnPa@4(Wy`Ww8<(A@;CC@kA
zDC#_)a`r0oZ|}K1&+csCT+XSsnr-%)EQ4w`19^?(Qd@4TY$z?1wo5y_ZJ}S}|EJ$h
z^jfU9@|;`TZv3dnI&!(%ianWPFUu3=f7<eC(WFqmQXYwUCf#p0y`HwA|H`U{JRJ|0
z8xNLPMg}gwn%$`+IN5+};?*4!!VGo@Fq!f_WV^t*<l=mF>3Jmvww7!KYt|Iyu%#}`
zxETMgGH=F;yin<X21!SL-COIwjc<;n&icljR~eDlYY%h;O3Nv+?yq6l62WTqE&GAC
zc2&GlVb39#vd+Md-p`~jhd$foc`&G{Bs{xVyyHjV$1Q&s+2*C!Hbtgyle#OjgC%UO
z*XqkJE?qH7ydv@Miu2nqKm5ETeurAG-28rn`<vr$*34D36wAJpY4DGK`N3%=5)K0S
z0n2R_CvJM(F!{<xEgr9_3mz$&2|Kg3o@@Q8rCjQhsk5!8$|s>FrQPY9^;$6lqn(`3
zFK>8tbLadpd*hgpt!aMpH$EsySLMWK-ktE~s4Q>FQ~5Np{HpzrIdU>BKDb+7IMV6o
z#XLpoP3d>zubUg2Ggh00-kzu6Z*Y5y_0hS0y$ecLWO?!ipFVr_N~ZKy{R}=0{#AS#
zCShm2&Q5Q+xU}#3`Ahm&jjMcGGXu<37Cu}db#-BH2D{8JtrY@$<9J^k;jmjG%W`(@
z@n0sftGcT!UTqL>mT&&3uwC`7NdL+R!F{{6&0j`e)KA)c+Ku&;TMzTpzQ@O_!nT~>
zH;sqIVFwRutSiU!)p~b1?&ouat^0HMkZ=E-O}xQ7?4O>}mOolI=kMEBd*@etx182j
z{bK3d)7jS_i*5_|e<-r%R`bbCwdZ#%J74km%+AOBCyJi*_a^LPw=8-Ts#|4#?edPS
z<+~LQ+l1}+Kk5<X>HlS4Pu@-?--_N#8CP%hzI;{UyK7tfvUW57?N=9^Gf11k7Ny-I
zx{S-acjK{;<j4gxmd7SUWW;4eF3dGrW4g|Ct(pD0TTAVf<jl6r+#g(<{Lk^W;=h!>
z-L`If<CLX(+(eUP(k5$Wgv?yH==uxZSqq|4!&Zl?Ee>29s4cPGx&O|h*(Ig&rNZS?
zL}V|_R#}|VH|gS}T_-doPgXn*R5o-gyx&|kr^_o=>sIQO-&ZC*xm<FY=d!_J(?Z^h
z3~t`xJ5*;ZzWOwp#c$DN#j-2Sf6H#&UU(x<e0^-uh574m%`Mv)tEXsPDwBD_C3gpd
zqS2GGpPQxrI_&dO7Zc;kpXOS!_>@xFj<bh7g!lJc`FMIq`uQ&D`R!*vp67cez9Q%|
zcS*tq{s=CfmD8Gz98U~bvtpkX&oraCacOq7g?%TpUQf(Fb?{4Yj6&(ko#kN<9cC*Y
zS($9yUoYkNwQYY+?^zAM*V{63@^W6~87iMOjh%lett4me?wpG@Mt`T5fB(8x!ZmrN
z&<5L=3j9$|9M*e>y2U-cz!ShDa!gCBWo4I?;|bAf$?C|@Z<Z#_+;}D`YyA_gsme#S
zlBbp~sZ}!%wv3#UIVW>&&VK&rDO`Tl=Ek$6*Y1ioZk~C_H}LF*6*+m+%y(SPDwgsO
zZPwJCal7|;hxkGL^dIlUV{W+Kmc4SpA#!?~io?f6N;hK4)lH2$<1fED9D05As=D8f
zyfRCwe%*|*?2=wO=Y(&_rX8yesr=v9JZb;88F`&cy~@_ry-v>A`~Q8m-jQ4Urt70_
z-r1BCYHaxM%h%6I^ULS_`OSQzdF%WCFaO%_|GHfGN?(J(lV<q?3%0fku!|Y)U;bxt
zns;8D()(++6PL3I+n2E&Jomy<fBW3s2R0s7*m=q@Ze3sYlSeY#s{&kSw3IC|J{fgs
z=H%(SWgjLl&^n>JYSP|BkF@rchgia${{1Z2_GU@o{a=nBs_Wj(zwo|G;7%;V-dYaL
z$MuuD9A~oM^nRoI^!c{+HFx57ytYdcu3=ny_i3<^&w<Bb5f8L7c8D{Sh0W!A)O%3Y
zQE1D}IlCUX{cJA2%2-{XyjLVf)l%RjC(FmomX&_WUxL~Kubb&g9Q_>qhdJ-aj?D)@
zGalcz>Tgu&@e7wURpYNWFVfOsxvTek?_0K?)y0lMCtsVS1)o@BxMM?Lwx`zj)KC8V
z&7E8m@_sG-DfZOQJ77W9znO2uv!2^Ni!FNea(Wcg)Y?~axw`XhDkAmGwp70O(fKs@
zH(%M3$;;Dv-fyi;Kc2j#^536_x8>eh^b6f%bnX?pq%e)A^$znr?dFmU&929R6-@K`
zr$6w0lKZJs|H<oTub<w2?B9F*<@=8r>~fu}7uFYUy3g{L)%^GW@9c-A=0C{&;rYj8
zkM*SYw?FLbzCZQ<{|DcW^BCqoN?h@iVc!7{2Kn$ti-g|)e=jmwFzB4z^=n()+4_Xn
z3~#ny31*nHJ@Ll``JEHRz1UJ%v}R7wzWJCj|IvT@{~sIX6);%r6rbq)%{K6lbi-+r
z{;NmyFW1eQ{GP#b77NS9V{tLf-xe)?VX=>E@wbJ}D@6Bl+XOVtbLqaYSm7a;313v$
z;)x%MF8tFdI}#lj{)MODgz*cN2lra1Ein7=?d`|>#>wFZ=?Pc39(06+@bI0~*}}))
zciwgC3#Q)+r`($_$_L2i^-Mp0{kZwj@<q=tDn7eg!CufK?A!D`<xlmy{Ez+t=@J_y
zb6v$hs6AyY$bIM3+wfClR_`{$P@O*i-v2z-9~S>k;Mb9A@_UplP{->Qe8K!ipSa=Q
zglW$oP5d#t;i+3_hSaULTLp4*4u4`&Ydrp_d<aZ8kx#qZHS3~d#o`yHUj(legg>@^
z9Nsv4^92P9f9+ok^2tfEr;}W#T<G{G5y@b?`GWjK>lY?`53Zlo@jIMaEVSLAPE2Km
zwG{L2T0ytE-g*`92eGI2oc-4SopoP)qujD9Zutw(Zxk*+diX=crOq=-{>Gaf&uvP`
zDL0Y4%CJ|y{igll#F_~iCp4c~eoOr(?&|d>DTdu#$jN`f@(ZscxYfTtWc$Nl*S=rn
z?}>^}EI*qnr`$bTcf!9R_<nQRqS6<sUzE)r&8uM9C+&Z@`-9dWgL_^6k1qc((o|b~
zvZSbW=CLcr>o?Xuls#&;fxk$5r)k&ejh}AJ-e`Zkt|$F?_#(jr{HGJX&%SuMV3YF;
z5i9=ptY284ZIE(Om3?!h|M-L74`Z9IKTQ9qw?}&WA+6;f#eYQq(f`BzueHHRk?WCC
z-6XF3hBIp4C%T`U{6u8W@{{%_{Cj#mQxEyb1)m5p*qji)dPeUDxdglPxo<Sh_)5+_
z+tx7g_R+m(WS_6oGWo;IBy{ip*PVT5`pze5d<)3&USoMp>c33MMCOXFtcgxH{<2hL
zKis#S;rLtUbcMRjhj~*T*;N<>q}i-H`gK9P{DISs(RTXp|BBYItoRaM!}@KFYF_2}
zSP!|E-fah~`DZ%4I_b8lbIYXSC+ThJS9<1cG?;mYVOiJALuV4j45DXne75_n`I)tN
z?oZKs7ORu?=4eSL9e%@5W@A0QcGlmvziq3Z{5w*+G(a;=D@<!qSN$j3A9jBV>R95B
z=0D*7l-`>4;GFpd#xJkd8ujZNot_c&j3LUuD1!ZGLk+_p8#&*!Hh~8ro4D;=#2$xD
zI(9bvi`w3g0VPv;8Q0!dyM40s(rn45(;qZKymzUm1<$xST|qnLzQfx$f)=;g4L;Yj
zvgs>^i%4DR5!7rjy`QveW8yRa7s}6=p508paC{Bd?E>q{iIRLVAEa)#YrEGy^ZddU
z+%LYEd5MV68T~_14<0Z4p?ZOPaiQ5ib76O$8~ekioUh@mJ6H0_{0d*uAAu?!|Fq?c
zzG|@RJKbL(?x1B=SRat0rSro!=&)Ji<||zsyQTfY8!j6ho_Vn-=4eiEjN<c)F(qB0
z3(wE6-S9j3iD()}^;GH3<&7*lZrXvhHIh*Wrf!yhw6f!1a^Qr$9Jg`>C!hVL7|#{8
zRH)>Z<Air+7nqNp5y;rc{G5kvp9Hg1lJuVb=6MN?%~=k|-ZSZbW6TQeQ42JgdO7ie
zl|aT{O*5h1?)?qhjxAqgf5d#@xdr@TN5v2J?7k45!F;PV{K&g@hO^9*K1<4-`KrEr
z$G?x~mUPZhUwmTP?f2p<WRH}6k=?U?!M_#!<r1bF>s=QI1lHvK(LZokazWzuMV&AF
zN;v&~+O9E}VA)ob<b3j9Mg20qbJdK!9s4D|CyCAyo6k7!vfvk+^9wm^BDD0Tb^W>%
zcw=|C<M}^(?#y}f?(VL#mvcnVy_qsmj%yif&|kat>PO!eS-lZ>Jj>UudByci)2yjM
z4@6~73-oOAV*b0jc49zj3ok1h-y5ITjhq(@9!_g})xkgO*Zn@>kD;eC_T4sAJ+SAy
zy}?qJWhKJG-TRzEZdjfUm-^u9H;ttx&?e8ta?16oELt@a%XZqh=;cp29i{R6ufc)v
ztodP&x*uf4XD<1_C)0Q7RaJw#bz4`(FSsuf`dnvA{L;2XQeJ}53-au@eLYyzc*Cpu
zt-8mBL(7!TIe2Cs5!6^{-C=O>66a*)!|qe2Puub5U)Q-v!8`du2}vxhzudkgO0s^P
zJiSxnmup3(WIV4eZ`9-4^Gnxc=7zWkuRc8MWaCL;t<%b9TDW$#De}CY=90s$f3tq|
zp(n@cp03V&TWaXL^qZ9YS~JtNpIpm&H1==3?ti(9YX<YNNeTa7E&I3OOO!&=O%3@q
zZebx)+&^$dTypn!@}FvPe~N=<?#gvvU8gkeYs~cBdbPmx+%gWSod@o}De-E5BtCKT
zmFp(EY##KB>|_!Z_PX=??261MO=9g&itpUJY?-}3!$Y8Kzu$bFj4j*6mfd+X+1toU
zcWUn*&q$wg5eXI^!H>S54*H#9E?73tWacGdQ{Pvd86QFtBt&QOm<I_5ep{U<RDH3`
z%I*;B(T4Re*;!-+GX5Wbt1CQ7w`;)<=0&cj)XgJ#f8|e6epsMlJb4qxeg1bM-)^|C
zdtQ2Ombhx8!`v6mJSW=Qm~+)bzP#dk^8NQq*7coU8A>5ay~?MAOXX+1$Z`w1Uf|ji
za(P~|d6iCbhX3Q<d9x3%W<GFoZ4pOeC&$9s^;H+MwB}yZZd&`l|D9&%zbEFh@tMJ;
zleVQU&R?zitNPv>?`cZQ>Yny=t9TzddSdT^DO1fSW-g5pPHg_Le8SACqAB7H+E2ya
zyp|VS&8~KdYwG>~^)e#8@k<v!d9yw7WUgCPYySH0`_|Uo{$}FKlOxfSaPvi>ukWhc
zzZ$Lzar1el9r&10vh)6R{ib6LDQpk<ou17T(TWk|^tjfq#PDBr)}*Vq1!7MYvD!?r
zwod=b?t1;o75fLT?}ywqQIyJ`Bk`fLahdZI3(I>U;Ze2Cg^NCJGhOHLf2*j&>dYPC
z=|BBW&i`NN^!3w=kG~GdvrJcr?nvM7bAHLsM(bJk|5nBy57{`~Kw_2Gf#aOQO#=7T
z%|yNmUp)Ui$uenu%I|#tRZO#PuzffueL>=rkXBm_hc2@P<1r?YrnmWvzsaA9-m!7T
zq-u6s-mj)b(oU`y>;#UjYp9ZX(!HQweM|ZUqb*A}$1jR)2xj}Ba6s#YG~a{vBF;K?
z+Z7A;IX;@nIJ>Q4VhY0*M~PkRVLwFwOxCgRRA%z$?^9P<D6^-#-eE(fnA4d&25W&!
zha&_FTr<oarX;`NERbTl-K!v?biky+UtrV0Ki^_FqS&2YaK7G=xK8ARZQ8H@UyP;(
zDQY%#v@kCZWqH3q|BcJyR}B?SBJ309ZCt3$P{_bmsrZ)5fG5G`09(Ub1}@3|ZD*Jw
zxC%rbth%FYxb2y*Wst}w)1qrZ*2>ZA{5(RmS3k)*9V^;dbE3l~&Pb#2ns(zgTl0%|
zmgXz@pZ9ulETK&KZ0M~w$w$MEb@?$1gs^L><ty19d+^{*f{q50gM{lQ*Q6fZ%!O~>
zx=V^Bwi!q;3FZHhXr8^b?M;$s=VWbxN!2I6FR{Jod@H)7?akYx<&*R;*flI>;Q7gP
z@BGQ<8Jaseu3T}BV%Z@3K;wXBgEj+Kk>EY<4|Q&K%xnAu!yA+uJ~G^9>|=^yd~kC4
zgsvyTRpv!rr(GFZYM54f*tI3s@~oU!!?Lvhg!_iA#%F)B-fVy3v?F`d170<5v+zv`
zThiQ5wx6{x{p5CXy6(OW<?e;%MPjZD96DuUYWCeV42Kz4F}`9D39w>cTD(FqLHPh*
z14jeXu{q6#&lApYHhgD1w=j=E=s*EO<D9&<)J0E1idU@PeuL}YZ)*Y1>XV+I*3H@P
ze{b4n`-$71{INY@J7K;;zIpYY6%o9DIgWKq=$=sCTGzMd-lViI$GUl+XiPZHoeA<%
zBEw@wGX{+lfs=08+{|wA<XbGP;2w9WXW!J%;ukhEP0HFEv_k&FjDEF+apex2liu!o
ztEkkvEvHq>*XgynW2f4}r`Zc#CSN#|8nDPS<Cu-cMHQ2sVqPYkz0&-uvs&7E+k{jC
zKYdx9`SRYM#LFAkNBF;E*1F3cvWtI3iF}fxMT_{~BbWBtUYs+NQ{=W<t95FX^zDMj
zf_E&0%<a0Rmp?eBx}#7sPv=PQ_Xo!e@=6}N@#jq`c$|1g*yCcy(rdR@u}+dTpLd|Q
zx-gkJkLTerg*>0by>@f-T|4(z2K3FDdMeV}aH{r+_b(=|ZMXf^xTr26+rHy(poc%>
zk*Bj*Joa5-*>u?X^Ob357P;*8UjHWYDDz2agAcqr+RQJoGKlHD;7wR;o-oNgfys7>
zE5l>PKE@cv4rUDoYf-a(j(%S!l(A|sN$_tFdBD}+bnl!fLsOdS_8tEBBCbDK)_zf{
zk~vPQ_IGEzm*Cx`2$>T(1<8Ba`9B}ARgsT;DQ>pd(ZFTnw+H>V>i@5BJ(Z*y>2J0`
zWJgEu$r|U+{+jnR>h@Tc+=$LO*ssXX9n;Poz`?HgzEkQ+m0S_G$$8hjPexB}o_)8`
zz~{cqsZZMyWo1h*{GR-K()Z}{DGnbu{oe3>{XUKE*?a37r%mUtI=?Ht`u=Mv(FIRE
zG`bYCKc8IDe#driwf~6~rxzqgZ!YIN%U317o4sQB?x%g4oqg3(3)EijRd#;7hdFVd
z(07;hR~U-rzKHF7|4~eB!>1QBn(yUpXf0SYgLC%VC)31vE1qy3f8zepM6>(;q3NFw
zD}<bXwz<Oo{nOwI@qDqKCXE*?y$gg`92GY(2~J2@_WzdLvST;a2xKlx$UJo_b>f97
zDmkVHnHnqfl9En&++4A5hfYzBkk;%M-7)f_9Z8}VCNT^2`y@^W+<dWb;=G5OH00$a
zPx4Q9{qeAT`upP*wX9p*PB;hh-187Qo_U!|oL81lv4>at+Rlu58<S#PSvN&L7roit
zQ_I8p*Sk)ArQfoL_vL=)KAm!6+NQv#6Q>7hDc#&6;J;F%+r}xm+xbWQkF3m-ihulP
zg-(%A`t_iGYTuH#KJzDab#2m+(|)Wd|4;n^<MBVi$AoQ8xPP4B|3T%$Hq{SJ267)R
zWw2FvsdAmFa9kR6EZOkr?fVnLqBhy<`X1Z(GoZQWa>D1kUByZq76%vfB<H%Qr_2$)
z=GG<rQhD)?yA6*DPix%hykpcQZJj!0<Nk#uzKhoE>)aF-*LiUN!yUpFx4O*>wum#T
z|EaHC_Jy~-(EELt=TsBD)qGR@7+3VInQ4AjE_2R{&u{04Oys%YbX8d0ck#~_&52zH
zrc64q#Gy04{{5-@FV`)(yO&w@cVOzv7oYx4y)<RgvwK1BFCA!^)HRLqaPjB-=<Kg&
zw-&YkKbLa&Wy;QywTwUaU!Lo~oBhJpPcNf>maJvG#89cL^7+SyPk)(SMKdNg%-Jrm
zmFvgWTTj=&$QFJ5>%^JcxAWs4Wihv>hzB``ENSG3=62X^C0?$#ELI_Si?@68`;_Um
zp85M%^S<{g>%5S|@x}F_WwG=#rRq#y&qwPS<uYa&PZzS8*e3mNJNqv$#Xh~v8O;Ha
zD#am@OE@0&{k7p=d-aRQ#5nU?S6&4b{1@(xU9eF+H)Y1`jxB2^Y}1_TdE|fP|KiCz
z|7hnbTi$6><)7Cqr8T2td5#HVsDYK_(=Cm=UCuf?v&w8>Y%e_AYxyL1v5Vxz&-WTE
zizhEtk$tA}`G|(ZrGjP|b9o8oRj>9lT3=AQ3|b!yU6Kl2ybD^>3R*F&prBxAU}0ea
zQUT&yVlIC*1g+R5?g-7;NoVt#4Fp==*Kx^Tc4m2aJWr)1WT`mc-@wau3^%;HLU@dx
z9Jqb`b!X7Z_9OD&&Q-YIdve9bJTD^k(>wLzWuf1iH_w~)@)hsNKD}n)cAr}Z*KA+^
zBmY>Ln5V$}?Ws#Q7XDRP<)tS7upsiUTDgDOs~3(Qnv>1*Px8N0Dvw`oyzSie%6|(d
z-q8vUynZQ2{!QGipIc=;)dKhA=UB*J%au>57oC6G$Ef>0`<4_n@6&grKCRt7Va25_
z(?aiGZ<_PBeA@E~C3}D7?OpN7S;P181{S}AnGepm@p<0}zY%^$bw`Tc0oG!^b2E(#
zj@hxuC$-fd``6kjBNjOI!n8HI4<a|D-%Q+}oXelf5t?}=d)@UH@x}k8zGv+VKkUBQ
z{~O0!<=@TRX?GVbdpb)wtFK0X)ATue9lkYf)=mkZ6WeUM*~#swq-a&}M)QiSZy_(@
z3;%Q76ipA`p{?^LVte?*^v6GUop0V*CGV!lmpF&9?(e~iwSLbt(?b9JWs?51R~8h1
zs4)zVHc*r*C@2_P7=RRjc!r=w*&qgJ6*y#BHE7MZOKNd)QD#9&ei3MuG<X@fUt(^m
zg1)Oyc%)aDwqsFdVvcVJT+%HuH#4VFK?5wV2@%doOfOb2HiB^S^Gci`3%+&D%uEz?
zjg3qcj0}t{6bucF&7ljwJxdaEGLs$i(sNQ348YnPi<3bM$1M#^!Ca85Kzv;zLqo8D
zb7Fx<YG!&y3D_90P;g0VZkU1vL_dh*0po!D<(8R~YNTKeUcqf=2MYwyf_BU>fEZVr
zR03Yi9#T}A3JM>vs8eEbD#$qO;gMNXT;iOOSOoLDPa>3WXoLt2kY6AntN@M!NRD#N
zOU_Tp%u83$56{eV%qz}>34pB)&rB)FC{~D8Ffs>u0hAG8*x1-u0fZrZ5DTORgiTG&
z6mVgPYHXyr5ojg6893EI*dSFPb3il<o0)-D{lj=58f2pZn)^%*Kr-lhjf{*GER2oe
zZZiP0KrRQV2bqJ=0TKaW@CE?|{gC|7yiCx6y$Yc6DToUZlbER)v<aaY9D3fFDWI?e
zxi|=#;+*qKLHiPr4FGu?mKL0gQWHz^ixf1H6ElnQ6%6$ZP4$cvG%`v`3as??!Gh^U
zi3J&%$;Ep4Md_N5gdb3ppHiBfN{Til4}dm3!0ZR5HwAs?#FE6E{B*D%z`lWQq(Eq@
zC`wJ^GEguy;WB^&1v67qV^f7R1&ElDk*SfT0$5fd4=QG0X=#ZeW?*1wimA@r3{A|)
z$iUnXP0Y~1(7*yi%ovoCKo-O8H8M0dK)1ul&=|vgMuz4V=4j>_8CrnSB&vBvrj{n?
z>Wsjt7N)MGC^0i9wTKH;$^~aurGml=v{NG}Kfgr55ESN+{O6gMmakw4N(10jRa}x-
UR00lNBO`NjQ!Z6iSARDy0N{{9KmY&$
literal 0
HcmV?d00001
diff --git a/example/random.c b/example/random.c
new file mode 100644
index 0000000..9e740c4
--- /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 0000000..fcff032
--- /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);
--
GitLab