mirror
/
libdspl-2.0
kopia lustrzana https://github.com/Dsplib/libdspl-2.0
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

fixed randgen bugs

pull/6/merge
Dsplib 2019-10-27 23:38:18 +03:00
rodzic 1e9430426d
commit 4e259f742f
9 zmienionych plików z 135 dodań i 82 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

3
dspl.workspace.layout
Wyświetl plik

@ -1,5 +1,6 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
<CodeBlocks_workspace_layout_file>
<FileVersion major="1" minor="0" />
<ActiveProject path="examples/examples.cbp" />
<ActiveProject path="lapack/lapack_double.cbp" />
<PreferredTarget name="Debug" />
</CodeBlocks_workspace_layout_file>

59
dspl/dspl.layout
Wyświetl plik

@ -2,31 +2,6 @@
<CodeBlocks_layout_file>
<FileVersion major="1" minor="0" />
<ActiveTarget name="Debug" />
<File name="..\include\dspl.c" open="1" top="0" tabpos="5" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="10057" topLine="192" />
</Cursor>
</File>
<File name="..\include\dspl.h" open="1" top="0" tabpos="4" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="8601" topLine="219" />
</Cursor>
</File>
<File name="src\fft.c" open="1" top="0" tabpos="3" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="1917" topLine="33" />
</Cursor>
</File>
<File name="src\randgen.c" open="1" top="0" tabpos="9" split="0" active="1" splitpos="0" zoom_1="5" zoom_2="0">
<Cursor>
<Cursor1 position="4569" topLine="53" />
</Cursor>
</File>
<File name="src\matrix.c" open="1" top="0" tabpos="6" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="4148" topLine="155" />
</Cursor>
</File>
<File name="src\inout.c" open="1" top="0" tabpos="7" split="0" active="1" splitpos="0" zoom_1="3" zoom_2="0">
<Cursor>
<Cursor1 position="1882" topLine="0" />
@ -37,9 +12,19 @@
<Cursor1 position="2457" topLine="53" />
</Cursor>
</File>
<File name="src\filter_ft.c" open="1" top="0" tabpos="8" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<File name="src\fft.c" open="1" top="0" tabpos="3" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="651" topLine="129" />
<Cursor1 position="1917" topLine="33" />
</Cursor>
</File>
<File name="..\include\dspl.c" open="1" top="0" tabpos="5" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="10057" topLine="192" />
</Cursor>
</File>
<File name="src\matrix.c" open="1" top="0" tabpos="6" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="4148" topLine="155" />
</Cursor>
</File>
<File name="src\cheby.c" open="1" top="0" tabpos="1" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
@ -47,4 +32,24 @@
<Cursor1 position="2438" topLine="60" />
</Cursor>
</File>
<File name="..\include\dspl.h" open="1" top="0" tabpos="4" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="48641" topLine="784" />
</Cursor>
</File>
<File name="src\mt19937.c" open="1" top="0" tabpos="13" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="5094" topLine="30" />
</Cursor>
</File>
<File name="src\filter_ft.c" open="1" top="0" tabpos="8" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="651" topLine="129" />
</Cursor>
</File>
<File name="src\randgen.c" open="1" top="0" tabpos="9" split="0" active="1" splitpos="0" zoom_1="3" zoom_2="0">
<Cursor>
<Cursor1 position="1170" topLine="26" />
</Cursor>
</File>
</CodeBlocks_layout_file>

49
dspl/src/mt19937.c
Wyświetl plik

@ -1,15 +1,15 @@
/*
/*
A C-program for MT19937-64 (2004/9/29 version).
Coded by Takuji Nishimura and Makoto Matsumoto.
This is a 64-bit version of Mersenne Twister pseudorandom number
generator.
Before using, initialize the state by using init_genrand64(seed)
Before using, initialize the state by using init_genrand64(seed)
or init_by_array64(init_key, key_length).
Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
@ -22,8 +22,8 @@
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
@ -40,12 +40,12 @@
References:
T. Nishimura, ``Tables of 64-bit Mersenne Twisters''
ACM Transactions on Modeling and
ACM Transactions on Modeling and
Computer Simulation 10. (2000) 348--357.
M. Matsumoto and T. Nishimura,
``Mersenne Twister: a 623-dimensionally equidistributed
uniform pseudorandom number generator''
ACM Transactions on Modeling and
ACM Transactions on Modeling and
Computer Simulation 8. (Jan. 1998) 3--30.
Any feedback is very welcome.
@ -68,25 +68,27 @@
/* initializes mt[NN] with a seed */
void mt19937_init_genrand64(unsigned long long seed, random_t* prnd)
{
unsigned long long *mt = prnd->mt19937_mt;
unsigned long long *mt = prnd->mt19937_mt;
int mti = prnd->mt19937_mti;
mt[0] = seed;
for (mti=1; mti<NN; mti++)
for (mti=1; mti<NN; mti++)
mt[mti] = (6364136223846793005ULL * (mt[mti-1] ^ (mt[mti-1] >> 62)) + mti);
prnd->mt19937_mti = mti;
}
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
void mt19937_init_by_array64(unsigned long long init_key[],
unsigned long long key_length,
unsigned long long key_length,
random_t* prnd)
{
unsigned long long i, j, k;
unsigned long long *mt = prnd->mt19937_mt;
unsigned long long *mt = prnd->mt19937_mt;
/* int mti = prnd->mt19937_mti; */
mt19937_init_genrand64(19650218ULL, prnd);
i=1; j=0;
k = (NN>key_length ? NN : key_length);
@ -100,15 +102,15 @@ void mt19937_init_by_array64(unsigned long long init_key[],
for (k=NN-1; k; k--) {
/* non linear */
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 2862933555777941757ULL)) - i;
i++;
if (i>=NN)
{
mt[0] = mt[NN-1];
i=1;
if (i>=NN)
{
mt[0] = mt[NN-1];
i=1;
}
}
mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */
mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */
}
@ -118,15 +120,15 @@ unsigned long long mt19937_genrand64_int64(random_t* prnd)
int i;
unsigned long long x;
static unsigned long long mag01[2]={0ULL, MATRIX_A};
unsigned long long *mt = prnd->mt19937_mt;
unsigned long long *mt = prnd->mt19937_mt;
int mti = prnd->mt19937_mti;
if (mti >= NN) { /* generate NN words at one time */
/* if init_genrand64() has not been called, */
/* a default initial seed is used */
if (mti == NN+1)
mt19937_init_genrand64(5489ULL, prnd);
if (mti == NN+1)
mt19937_init_genrand64(5489ULL, prnd);
for (i=0;i<NN-MM;i++) {
x = (mt[i]&UM)|(mt[i+1]&LM);
@ -149,6 +151,7 @@ unsigned long long mt19937_genrand64_int64(random_t* prnd)
x ^= (x << 37) & 0xFFF7EEE000000000ULL;
x ^= (x >> 43);
prnd->mt19937_mti = mti;
return x;
}

21
dspl/src/randgen.c
Wyświetl plik

@ -32,24 +32,35 @@
/******************************************************************************
random generator initialization
*******************************************************************************/
int DSPL_API random_init(random_t* prnd, int type)
int DSPL_API random_init(random_t* prnd, int type, void* seed)
{
srand(time(NULL));
if(!prnd)
return RES_OK;
switch(type)
{
case RAND_TYPE_MRG32K3A:
/* MRG32k3a init */
prnd->mrg32k3a_x[0] = prnd->mrg32k3a_x[1] = 1.0;
prnd->mrg32k3a_y[0] = prnd->mrg32k3a_y[1] = prnd->mrg32k3a_y[2] = 1.0;
prnd->mrg32k3a_x[2] = rand();
if(seed)
prnd->mrg32k3a_x[2] = *((double*)seed);
else
prnd->mrg32k3a_x[2] = (double) rand() * rand();
break;
case RAND_TYPE_MT19937:
mt19937_init_genrand64((unsigned long long)rand()*rand(), prnd);
if(seed)
mt19937_init_genrand64(*((unsigned long long*)seed), prnd);
else
mt19937_init_genrand64((unsigned long long)rand()*rand(), prnd);
break;
default:
return ERROR_RAND_TYPE;
}
prnd->type = type;
return RES_OK;
}
@ -118,7 +129,7 @@ int DSPL_API randu(double* x, int n, random_t* prnd)
{
int i;
if(!x || !prnd)
if(!x)
return ERROR_PTR;
if(n < 0)
return ERROR_SIZE;
@ -144,7 +155,7 @@ int DSPL_API randu(double* x, int n, random_t* prnd)
if(n<1)
return ERROR_SIZE;
for(i = 0; i < n; i++)
x[n] = (double)rand()/RAND_MAX;
x[i] = (double)rand()/RAND_MAX;
}
return RES_OK;

15
examples/examples.depend
Wyświetl plik

@ -20,3 +20,18 @@
<string.h>
"dspl.h"
1572196728 source:f:\dsplib.org\libdspl-2.0\include\dspl.c
<windows.h>
<dlfcn.h>
<stdio.h>
"dspl.h"
1572200729 f:\dsplib.org\libdspl-2.0\include\dspl.h
<math.h>
1572200804 source:f:\dsplib.org\libdspl-2.0\examples\src\filter_iir_test.c
<stdio.h>
<stdlib.h>
<string.h>
"dspl.h"

2
examples/examples.layout
Wyświetl plik

@ -14,7 +14,7 @@
</File>
<File name="src\filter_iir_test.c" open="1" top="1" tabpos="1" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="936" topLine="8" />
<Cursor1 position="826" topLine="28" />
</Cursor>
</File>
</CodeBlocks_layout_file>

36
examples/src/filter_iir_test.c
Wyświetl plik

@ -4,10 +4,14 @@
#include "dspl.h"
#define ORD 6
#define N 1000
#define N 2000
int main(int argc, char* argv[])
{
void* handle; /* DSPL handle */
void* handle; /* DSPL handle */
double b[ORD+1], a[ORD+1];
double t[N], s[N], n[N], sf[N];
@ -15,30 +19,38 @@ int main(int argc, char* argv[])
int k;
int err;
handle = dspl_load(); /* Load DSPL function */
/* Load DSPL function */
handle = dspl_load();
random_init(&rnd, RAND_TYPE_MRG32K3A); /* random generator init */
linspace(0, N, N, DSPL_PERIODIC, t); /* fill t vector */
randn(n, N, 0, 1.0, &rnd); /* generate noise */
/* random generator init */
random_init(&rnd, RAND_TYPE_MT19937, NULL);
/* fill time vector */
linspace(0, N, N, DSPL_PERIODIC, t);
/* generate noise */
randn(n, N, 0, 1.0, &rnd);
/* input signal s = sin(2*pi*t) + n(t) */
for(k = 0; k < N; k++)
s[k] = sin(M_2PI*0.02*t[k]) + n[k];
/* iir filter calculation */
/* IIR filter coefficients calculation */
iir(1.0, 70.0, ORD, 0.06, 0.0, DSPL_FILTER_ELLIP | DSPL_FILTER_LPF, b, a);
// input signal filtration
/* input signal filtration */
filter_iir(b, a, ORD, s, N, sf);
// save input signal and filter output to the txt-files
/* save input signal and filter output to the txt-files */
writetxt(t,s, N, "dat/s.txt");
writetxt(t,sf,N, "dat/sf.txt");
/* run GNUPLOT script */
err = gnuplot_script(argc, argv, "gnuplot/filter_iir.plt");
dspl_free(handle); // free dspl handle
// run GNUPLOT script
/* free DSPL handle */
dspl_free(handle);
return err;
}

27
include/dspl.h
Wyświetl plik

@ -58,17 +58,17 @@ typedef struct
#define RAND_MT19937_NN 312
typedef struct
{
double mrg32k3a_seed;
double mrg32k3a_x[3];
double mrg32k3a_y[3];
/* The array for the MT19937 state vector */
unsigned long long mt19937_mt[RAND_MT19937_NN];
unsigned long long mt19937_mt[RAND_MT19937_NN];
int mt19937_mti;
int type;
}random_t;
@ -495,7 +495,7 @@ DECLARE_FUNC(void, fft_free, fft_t*);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, fft_mag, double* x
COMMA int n
COMMA fft_t* pfft
COMMA fft_t* pfft
COMMA double fs
COMMA int flag
COMMA double* mag
@ -503,7 +503,7 @@ DECLARE_FUNC(int, fft_mag, double* x
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, fft_mag_cmplx, complex_t* x
COMMA int n
COMMA fft_t* pfft
COMMA fft_t* pfft
COMMA double fs
COMMA int flag
COMMA double* mag
@ -722,7 +722,7 @@ DECLARE_FUNC(int, low2low, double* b
COMMA double* beta
COMMA double* alpha);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, matrix_eig_cmplx, complex_t* a
DECLARE_FUNC(int, matrix_eig_cmplx, complex_t* a
COMMA int n
COMMA complex_t* v
COMMA int* info);
@ -745,13 +745,13 @@ DECLARE_FUNC(int, matrix_mul, double* a
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, matrix_print, double* a
COMMA int n
COMMA int m
COMMA int m
COMMA const char* name
COMMA const char* format);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, matrix_print_cmplx, complex_t* a
COMMA int n
COMMA int m
COMMA int n
COMMA int m
COMMA const char* name
COMMA const char* format);
/*----------------------------------------------------------------------------*/
@ -799,7 +799,8 @@ DECLARE_FUNC(int, randn, double*
COMMA random_t* prnd);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, random_init, random_t* prnd
COMMA int type);
COMMA int type
COMMA void* seed);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, randu, double*
COMMA int
@ -879,7 +880,7 @@ DECLARE_FUNC(int, verif, double* x
COMMA double* y
COMMA size_t n
COMMA double eps
COMMA double* err);
COMMA double* err);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, verif_cmplx, complex_t* x
COMMA complex_t* y

5
lapack/lapack_double.layout
Wyświetl plik

@ -2,4 +2,9 @@
<CodeBlocks_layout_file>
<FileVersion major="1" minor="0" />
<ActiveTarget name="Debug" />
<File name="src\dgbbrd.f" open="0" top="0" tabpos="0" split="0" active="1" splitpos="0" zoom_1="0" zoom_2="0">
<Cursor>
<Cursor1 position="16592" topLine="517" />
</Cursor>
</File>
</CodeBlocks_layout_file>