diff --git a/Makefile b/Makefile
index 3b58ffe..5b70a4e 100644
--- a/Makefile
+++ b/Makefile
@@ -6,8 +6,11 @@ all:
$(MAKE) -C lapack
$(MAKE) -C dspl
$(MAKE) -C examples
-
+ $(MAKE) -C performance
+ $(MAKE) -C verification
clean:
$(MAKE) -C dspl clean
- $(MAKE) -C examples clean
\ No newline at end of file
+ $(MAKE) -C examples clean
+ $(MAKE) -C performance clean
+ $(MAKE) -C verification clean
\ No newline at end of file
diff --git a/dspl/src/inout.c b/dspl/src/inout.c
index 6e4c215..b3cc5be 100644
--- a/dspl/src/inout.c
+++ b/dspl/src/inout.c
@@ -44,9 +44,6 @@ void DSPL_API dspl_info()
}
-
-
-
#ifdef DOXYGEN_ENGLISH
#endif
@@ -90,7 +87,6 @@ int DSPL_API readbin(char* fn, void** x, int* k, int* dtype)
if(k)
*k = n*m;
-
switch(t)
{
case DAT_DOUBLE:
@@ -111,6 +107,7 @@ int DSPL_API readbin(char* fn, void** x, int* k, int* dtype)
err = ERROR_FREAD_SIZE;
goto exit_label;
}
+ break;
default:
err = ERROR_DAT_TYPE;
goto exit_label;
@@ -323,6 +320,7 @@ int DSPL_API writebin(void* x, int n, int dtype, char* fn)
}
+
if(fwrite(&n, sizeof(int), 1, pFile) != 1)
{
res = ERROR_FWRITE_SIZE;
diff --git a/dspl/src/statistic.c b/dspl/src/statistic.c
index 4d5cf46..05991cb 100644
--- a/dspl/src/statistic.c
+++ b/dspl/src/statistic.c
@@ -257,4 +257,8 @@ int DSPL_API minmax(double* x, int n, double* xmin, double* xmax)
*xmax = max;
return RES_OK;
-}
\ No newline at end of file
+}
+
+
+
+
diff --git a/dspl/src/xcorr.c b/dspl/src/xcorr.c
new file mode 100644
index 0000000..0a9a9eb
--- /dev/null
+++ b/dspl/src/xcorr.c
@@ -0,0 +1,308 @@
+/*
+* Copyright (c) 2015-2019 Sergey Bakhurin
+* Digital Signal Processing Library [http://dsplib.org]
+*
+* This file is part of DSPL.
+*
+* is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* DSPL 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 General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with Foobar. If not, see .
+*/
+
+#include
+#include
+#include
+#include "dspl.h"
+
+
+int xcorr_fft_size(int nx, int ny, int* pnfft, int* pndata);
+
+int xcorr_get_lag_cmplx(complex_t* x, int nd, int nr, complex_t* r, double* t);
+
+int xcorr_krn(complex_t* x, int nx, complex_t* y, int ny, fft_t* pfft,
+ int flag, int nr, complex_t* r, double* t);
+
+int xcorr_scale_cmplx(complex_t* x, int nd, int flag);
+
+
+
+#ifdef DOXYGEN_ENGLISH
+
+#endif
+#ifdef DOXYGEN_RUSSIAN
+
+#endif
+int DSPL_API xcorr(double* x, int nx, double* y, int ny,
+ int flag, int nr, double* r, double* t)
+{
+ fft_t fft = {0};
+ int err;
+ complex_t *cr = (complex_t*)malloc((2 * nr + 1) * sizeof(complex_t));
+ complex_t *cx = (complex_t*)malloc( nx * sizeof(complex_t));
+ complex_t *cy = (complex_t*)malloc( ny * sizeof(complex_t));
+
+ err = re2cmplx(x, nx, cx);
+ if(err != RES_OK)
+ goto exit_label;
+ err = re2cmplx(y, ny, cy);
+ if(err != RES_OK)
+ goto exit_label;
+ err = xcorr_krn(cx, nx, cy, ny, &fft, flag, nr, cr, t);
+ if(err != RES_OK)
+ goto exit_label;
+
+ err = cmplx2re(cr, 2*nr+1, r, NULL);
+
+exit_label:
+ if(cr)
+ free(cr);
+ fft_free(&fft);
+ return err;
+}
+
+
+
+
+
+#ifdef DOXYGEN_ENGLISH
+
+#endif
+#ifdef DOXYGEN_RUSSIAN
+
+#endif
+int DSPL_API xcorr_cmplx(complex_t* x, int nx, complex_t* y, int ny,
+ int flag, int nr, complex_t* r, double* t)
+{
+ fft_t fft = {0};
+ int err;
+ err = xcorr_krn(x, nx, y, ny, &fft, flag, nr, r, t);
+ fft_free(&fft);
+ return err;
+}
+
+
+
+
+
+int xcorr_get_lag_cmplx(complex_t* x, int nd, int nr, complex_t* r, double* t)
+{
+ int i;
+ if(!x || !r)
+ return ERROR_PTR;
+ if(nd < 1 || nr < 1)
+ return ERROR_SIZE;
+
+ if(nr < nd)
+ memcpy(r, x+nd-1-nr, (2*nr+1)*sizeof(complex_t));
+ else
+ {
+ memset(r, 0, (2*nr+1) * sizeof(complex_t));
+ memcpy(r + nr - nd + 1, x, (2*nd-1)*sizeof(complex_t));
+ }
+ if(t)
+ for(i = 0; i < 2*nr+1; i++)
+ t[i] = (double)i - (double)nr;
+ return RES_OK;
+}
+
+
+
+
+
+#ifdef DOXYGEN_ENGLISH
+
+#endif
+#ifdef DOXYGEN_RUSSIAN
+
+#endif
+int xcorr_krn(complex_t* x, int nx, complex_t* y, int ny, fft_t* pfft,
+ int flag, int nr, complex_t* r, double* t)
+{
+ complex_t *px = NULL;
+ complex_t *py = NULL;
+ complex_t *pc = NULL;
+ complex_t *pX = NULL;
+ complex_t *pY = NULL;
+ complex_t *pC = NULL;
+
+ int nfft, ndata;
+ int err, i;
+
+ if(!x || !y || !r)
+ return ERROR_PTR;
+ if(nx < 1 || ny < 1 || nr < 1)
+ return ERROR_SIZE;
+
+ err = xcorr_fft_size(nx, ny, &nfft, &ndata);
+ if(err!= RES_OK)
+ goto exit_label;
+
+ /* memory allocation */
+ px = (complex_t*)malloc(nfft * sizeof(complex_t));
+ py = (complex_t*)malloc(nfft * sizeof(complex_t));
+ pc = (complex_t*)malloc(nfft * sizeof(complex_t));
+
+ pX = (complex_t*)malloc(nfft * sizeof(complex_t));
+ pY = (complex_t*)malloc(nfft * sizeof(complex_t));
+ pC = (complex_t*)malloc(nfft * sizeof(complex_t));
+
+ memset(px, 0, nfft * sizeof(complex_t));
+ memset(py, 0, nfft * sizeof(complex_t));
+
+ memcpy(px + ndata - 1, x, nx * sizeof(complex_t));
+ memcpy(py, y, ny * sizeof(complex_t));
+
+ err = fft_cmplx(px, nfft, pfft, pX);
+ if(err!= RES_OK)
+ goto exit_label;
+
+ err = fft_cmplx(py, nfft, pfft, pY);
+ if(err!= RES_OK)
+ goto exit_label;
+
+ for(i = 0; i < nfft; i++)
+ {
+ RE(pC[i]) = CMCONJRE(pX[i], pY[i]);
+ IM(pC[i]) = CMCONJIM(pX[i], pY[i]);
+ }
+
+ err = ifft_cmplx(pC, nfft, pfft, pc);
+ if(err!= RES_OK)
+ goto exit_label;
+
+ err = xcorr_scale_cmplx(pc, ndata, flag);
+ if(err!= RES_OK)
+ goto exit_label;
+
+ err = xcorr_get_lag_cmplx(pc, ndata, nr, r, t);
+
+exit_label:
+ if(px)
+ free(px);
+ if(py)
+ free(py);
+ if(pc)
+ free(pc);
+ if(pX)
+ free(pX);
+ if(pY)
+ free(pY);
+ if(pC)
+ free(pC);
+ return err;
+}
+
+
+
+
+
+/*
+Возвращает размер FFT для расчета полного вектора автокорреляции
+или кросскорреляции.
+
+Размер кросскорреляции равен
+N = 2 * nx - 1, если nx > ny;
+N = 2 * ny - 1, eсли nx <= ny.
+
+Посколку N может оказаться неудачным размером для FFT, то можно добить нулями
+до удобной длины.
+
+Если например N = 1025, то добивать до длины 2048 не очень эффективно, потому
+что много лишних нулей.
+
+Если мы рассмотрим N = 2^L + D, то целесообразно использовать
+
+ NFFT = 2^L + 2^(L - P), где P = 0,1,2 или 3.
+
+Тогда NFFT = 2^(L-P) * (2^P + 1). Тогда 2^(L-P) реалиуем как radix-2, а
+дополнительный составной множитель при P = 0,1,2 или 3 равен соответсвенно
+9, 5, 3 или 2, а для этих длин существуют хорошие процедуры.
+При P = 4 составной множитель будет (2^P + 1) = 17, что не очень хорошо.
+*/
+int xcorr_fft_size(int nx, int ny, int* pnfft, int* pndata)
+{
+ int nfft, nfft2, r2, dnfft;
+
+ if(nx < 1 || ny < 1)
+ return ERROR_SIZE;
+ if(!pnfft || !pndata)
+ return ERROR_PTR;
+
+
+
+ if(nx > ny)
+ {
+ nfft = 2*nx - 1;
+ *pndata = nx;
+ }
+ else
+ {
+ nfft = 2*ny - 1;
+ *pndata = ny;
+ }
+ nfft2 = nfft;
+
+ r2 = 0;
+ while(nfft2 >>= 1)
+ r2++;
+
+ if(r2 > 3)
+ {
+ dnfft = 1 << (r2 - 3);
+ while(((1 << r2) + dnfft) < nfft)
+ dnfft <<= 1;
+ nfft = (1 << r2) + dnfft;
+ }
+
+ *pnfft = nfft;
+ return RES_OK;
+}
+
+
+
+
+
+
+int xcorr_scale_cmplx(complex_t* x, int nd, int flag)
+{
+ int i;
+ double w;
+ if(!x)
+ return ERROR_PTR;
+ if(nd < 1)
+ return ERROR_SIZE;
+
+ switch(flag)
+ {
+ case DSPL_XCORR_NOSCALE:
+ break;
+ case DSPL_XCORR_BIASED:
+ for(i = 0; i < 2 * nd - 1; i++)
+ {
+ w = 1.0 / (double)nd;
+ RE(x[i]) *= w;
+ IM(x[i]) *= w;
+ }
+ break;
+ case DSPL_XCORR_UNBIASED:
+ for(i = 1; i < 2 * nd - 1; i++)
+ {
+ w = 1.0 / ((double)nd - fabs((double)(i - nd)));
+ RE(x[i-1]) *= w;
+ IM(x[i-1]) *= w;
+ }
+ break;
+ default:
+ return ERROR_XCORR_FLAG;
+ }
+ return RES_OK;
+}
diff --git a/examples/src/xcorr_cmplx_test.c b/examples/src/xcorr_cmplx_test.c
new file mode 100644
index 0000000..47f0002
--- /dev/null
+++ b/examples/src/xcorr_cmplx_test.c
@@ -0,0 +1,78 @@
+#include
+#include
+#include
+#include "dspl.h"
+
+#define N 200
+#define P 20
+
+int main(int argc, char* argv[])
+{
+ void* hdspl; /* DSPL handle */
+ void* hplot; /* GNUPLOT handles */
+
+ random_t rnd = {0}; /* random structure */
+ fft_t fft = {0};
+
+ hdspl = dspl_load(); /* Load DSPL function */
+
+ complex_t x[N];
+ complex_t z[2*P+1];
+ double t[2*P+1];
+
+ int err;
+
+ err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
+ printf("err = %.8x\n", err);
+ if(err != RES_OK)
+ goto exit_label;
+
+ err = randn((double*)x, 2*N, 0.0, 1.0, &rnd);
+ printf("err = %.8x\n", err);
+ if(err != RES_OK)
+ goto exit_label;
+
+ err = xcorr_cmplx(x, N, x, N, DSPL_XCORR_BIASED, P, z, t);
+ printf("err = %.8x\n", err);
+ writetxt_cmplx_re(t, z, 2*P+1, "dat/xcorr_b_re.txt");
+ writetxt_cmplx_im(t, z, 2*P+1, "dat/xcorr_b_im.txt");
+
+ err = xcorr_cmplx(x, N, x, N, DSPL_XCORR_UNBIASED, P, z, t);
+ printf("err = %.8x\n", err);
+ writetxt_cmplx_re(t, z, 2*P+1, "dat/xcorr_u_re.txt");
+ writetxt_cmplx_im(t, z, 2*P+1, "dat/xcorr_u_im.txt");
+
+ /**************************************************************************/
+ /* plotting by GNUPLOT */
+ /**************************************************************************/
+ /* Create window plot */
+ err = gnuplot_create(argc, argv, 420, 420, "img/xcorr_re.png", &hplot);
+ if(err != RES_OK)
+ goto exit_label;
+
+ gnuplot_cmd(hplot, "set xlabel 't'");
+ gnuplot_cmd(hplot, "set ylabel 'z'");
+ gnuplot_cmd(hplot, "unset key");
+ gnuplot_cmd(hplot, "set title 'xcorr'");
+ gnuplot_cmd(hplot, "plot 'dat/xcorr_u_re.txt' with lines,\\");
+ gnuplot_cmd(hplot, " 'dat/xcorr_b_re.txt' with lines");
+ gnuplot_close(hplot);
+
+ err = gnuplot_create(argc, argv, 420, 420, "img/xcorr_im.png", &hplot);
+ if(err != RES_OK)
+ goto exit_label;
+
+ gnuplot_cmd(hplot, "set xlabel 't'");
+ gnuplot_cmd(hplot, "set ylabel 'z'");
+ gnuplot_cmd(hplot, "unset key");
+ gnuplot_cmd(hplot, "set title 'xcorr'");
+ gnuplot_cmd(hplot, "plot 'dat/xcorr_u_im.txt' with lines,\\");
+ gnuplot_cmd(hplot, " 'dat/xcorr_b_im.txt' with lines");
+ gnuplot_close(hplot);
+
+exit_label:
+ /* free dspl handle */
+ dspl_free(hdspl);
+ return err;
+}
+
diff --git a/examples/src/xcorr_test.c b/examples/src/xcorr_test.c
new file mode 100644
index 0000000..c9e1654
--- /dev/null
+++ b/examples/src/xcorr_test.c
@@ -0,0 +1,63 @@
+#include
+#include
+#include
+#include "dspl.h"
+
+#define N 500
+#define P 490
+
+int main(int argc, char* argv[])
+{
+ void* hdspl; /* DSPL handle */
+ void* hplot; /* GNUPLOT handles */
+ random_t rnd = {0}; /* random structure */
+ hdspl = dspl_load(); /* Load DSPL function */
+
+ double x[N];
+ double z[2*P+1];
+ double t[2*P+1];
+
+ int err;
+
+ err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
+ if(err != RES_OK)
+ goto exit_label;
+
+ err = randn(x, N, 0.0, 1.0, &rnd);
+ if(err != RES_OK)
+ goto exit_label;
+
+ err = xcorr(x, N, x, N, DSPL_XCORR_UNBIASED, P, z, t);
+
+ /* Save to files "dat/randu_mrg32k3a.txt" */
+ writetxt(t, z, 2*P+1, "dat/xcorr_unbiased.txt");
+
+ err = xcorr(x, N, x, N, DSPL_XCORR_BIASED, P, z, t);
+
+ /* Save to files "dat/randu_mrg32k3a.txt" */
+ writetxt(t, z, 2*P+1, "dat/xcorr_biased.txt");
+
+ /**************************************************************************/
+ /* plotting by GNUPLOT */
+ /**************************************************************************/
+ /* Create window plot */
+ err = gnuplot_create(argc, argv, 420, 420, "img/xcorr_test.png", &hplot);
+ if(err != RES_OK)
+ goto exit_label;
+
+ gnuplot_cmd(hplot, "set xlabel 't'");
+ gnuplot_cmd(hplot, "set ylabel 'z'");
+ gnuplot_cmd(hplot, "unset key");
+ gnuplot_cmd(hplot, "set title 'xcorr'");
+ gnuplot_cmd(hplot, "plot 'dat/xcorr_unbiased.txt' with lines, 'dat/xcorr_biased.txt' with lines");
+
+
+exit_label:
+ if(hplot)
+ gnuplot_close(hplot);
+
+ /* free dspl handle */
+ dspl_free(hdspl);
+ return 0;
+}
+
diff --git a/include/dspl.c b/include/dspl.c
index 925774f..af3a482 100644
--- a/include/dspl.c
+++ b/include/dspl.c
@@ -185,6 +185,9 @@ p_writetxt_cmplx_im writetxt_cmplx_im ;
p_writetxt_cmplx_re writetxt_cmplx_re ;
p_writetxt_int writetxt_int ;
+p_xcorr xcorr ;
+p_xcorr_cmplx xcorr_cmplx ;
+
#ifdef WIN_OS
#define LOAD_FUNC(fn) \
@@ -488,6 +491,10 @@ void* dspl_load()
LOAD_FUNC(writetxt_cmplx_re);
LOAD_FUNC(writetxt_int);
+ LOAD_FUNC(xcorr);
+ LOAD_FUNC(xcorr_cmplx);
+
+
#ifdef WIN_OS
return (void*)handle;
exit_label:
diff --git a/include/dspl.h b/include/dspl.h
index 29b99e5..a16f792 100644
--- a/include/dspl.h
+++ b/include/dspl.h
@@ -561,6 +561,7 @@ Variable `z = 1-2j`, here `j` - imaginary unit, but variable `y = 5`.
#define ERROR_WIN_SYM 0x23091925
#define ERROR_WIN_TYPE 0x23092025
/* X 0x24xxxxxx*/
+#define ERROR_XCORR_FLAG 0x24031518
/* Y 0x25xxxxxx*/
/* Z 0x26xxxxxx*/
@@ -619,6 +620,11 @@ Variable `z = 1-2j`, here `j` - imaginary unit, but variable `y = 5`.
#define DSPL_FILTER_ELLIP 0x00000800
+#define DSPL_XCORR_NOSCALE 0x00000000
+#define DSPL_XCORR_BIASED 0x00000001
+#define DSPL_XCORR_UNBIASED 0x00000002
+
+
#define ELLIP_ITER 16
#define ELLIP_MAX_ORD 24
@@ -1400,6 +1406,24 @@ DECLARE_FUNC(int, writetxt_int, int*
COMMA int
COMMA char*);
/*----------------------------------------------------------------------------*/
+DECLARE_FUNC(int, xcorr, double* x
+ COMMA int nx
+ COMMA double* y
+ COMMA int ny
+ COMMA int flag
+ COMMA int nr
+ COMMA double* r
+ COMMA double* t);
+/*----------------------------------------------------------------------------*/
+DECLARE_FUNC(int, xcorr_cmplx, complex_t* x
+ COMMA int nx
+ COMMA complex_t* y
+ COMMA int ny
+ COMMA int flag
+ COMMA int nr
+ COMMA complex_t* r
+ COMMA double* t);
+/*----------------------------------------------------------------------------*/
#ifdef __cplusplus
diff --git a/performance/Makefile b/performance/Makefile
new file mode 100644
index 0000000..8356fc2
--- /dev/null
+++ b/performance/Makefile
@@ -0,0 +1,43 @@
+
+include ../make.inc
+
+
+SRC_DIR = src
+OBJ_DIR = obj
+
+# C-compiler flags
+CFLAGS = -c -O3 -I$(INC_DIR) -D$(DEF_OS)
+
+
+DSPL_C = ../include/dspl.c
+DSPL_O = obj/dspl.o
+
+
+# DSPL src and obj files list
+SRC_FILES = $(wildcard $(SRC_DIR)/*.c)
+OBJ_FILES = $(addprefix $(OBJ_DIR)/,$(notdir $(SRC_FILES:.c=.o)))
+EXE_FILES = $(addprefix bin/,$(notdir $(SRC_FILES:.c=.exe)))
+
+
+all: $(EXE_FILES) bin/$(LIB_NAME)
+
+
+
+
+bin/%.exe:$(OBJ_DIR)/%.o $(DSPL_O)
+ $(CC) $< $(DSPL_O) -o $@ $(LFLAGS)
+
+
+$(OBJ_DIR)/%.o:$(SRC_DIR)/%.c
+ $(CC) $(CFLAGS) $< -o $@ $(LFLAGS)
+
+$(DSPL_O):$(DSPL_C)
+ $(CC) $(CFLAGS) $(DSPL_C) -o $(DSPL_O)
+
+bin/$(LIB_NAME):
+ cp $(RELEASE_DIR)/$(LIB_NAME) bin/$(LIB_NAME)
+
+clean:
+ rm -f $(OBJ_DIR)/*.o
+ rm -f bin/*.exe
+
diff --git a/performance/bin/.gitignore b/performance/bin/.gitignore
new file mode 100644
index 0000000..a5a6f96
--- /dev/null
+++ b/performance/bin/.gitignore
@@ -0,0 +1,2 @@
+*.exe
+*.dll
\ No newline at end of file
diff --git a/performance/obj/.gitignore b/performance/obj/.gitignore
new file mode 100644
index 0000000..444cfd9
--- /dev/null
+++ b/performance/obj/.gitignore
@@ -0,0 +1 @@
+.*o
\ No newline at end of file
diff --git a/performance/src/xcorr_cmplx_performance.c b/performance/src/xcorr_cmplx_performance.c
new file mode 100644
index 0000000..1a7d2ae
--- /dev/null
+++ b/performance/src/xcorr_cmplx_performance.c
@@ -0,0 +1,57 @@
+#include
+#include
+#include
+#include
+#include "dspl.h"
+
+#define N 16777216
+#define P 16777216
+
+int main(int argc, char* argv[])
+{
+ void* hdspl; /* DSPL handle */
+ void* hplot; /* GNUPLOT handles */
+
+ random_t rnd = {0}; /* random structure */
+
+ hdspl = dspl_load(); /* Load DSPL function */
+
+ complex_t *x = (complex_t*) malloc (N * sizeof(complex_t));
+ complex_t *z = (complex_t*) malloc ((2*P+1) * sizeof(complex_t));
+ complex_t *y = (complex_t*) malloc ((2*P+1) * sizeof(complex_t));
+
+ int err, i, j, k = 1, size = N;
+ clock_t t;
+ double dt;
+
+ err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
+ if(err != RES_OK)
+ goto exit_label;
+
+ err = randn((double*)x, 2*N, 0.0, 1.0, &rnd);
+ if(err != RES_OK)
+ goto exit_label;
+
+ for(j = 0; j < 22; j++)
+ {
+ t = clock();
+ for(i = 0; i < k; i++)
+ xcorr_cmplx(x, size, x, size, DSPL_XCORR_BIASED, P, y, NULL);
+ t = clock() - t;
+ dt = ((double)t)/CLOCKS_PER_SEC/(double)k;
+ printf ("N = %8d xcorr time: %12.6f sec\n", size, dt);
+
+ k *= 1.7;
+ size /= 2;;
+ }
+
+
+exit_label:
+ free(x);
+ free(y);
+ free(z);
+ /* free dspl handle */
+ dspl_free(hdspl);
+ return err;
+}
+
diff --git a/verification/Makefile b/verification/Makefile
new file mode 100644
index 0000000..28513b2
--- /dev/null
+++ b/verification/Makefile
@@ -0,0 +1,42 @@
+
+include ../make.inc
+
+
+SRC_DIR = src
+OBJ_DIR = obj
+
+# C-compiler flags
+CFLAGS = -c -O3 -I$(INC_DIR) -D$(DEF_OS)
+
+
+DSPL_C = ../include/dspl.c
+DSPL_O = obj/dspl.o
+
+
+# DSPL src and obj files list
+SRC_FILES = $(wildcard $(SRC_DIR)/*.c)
+OBJ_FILES = $(addprefix $(OBJ_DIR)/,$(notdir $(SRC_FILES:.c=.o)))
+EXE_FILES = $(addprefix bin/,$(notdir $(SRC_FILES:.c=.exe)))
+
+
+all: $(EXE_FILES) bin/$(LIB_NAME)
+
+
+
+
+bin/%.exe:$(OBJ_DIR)/%.o $(DSPL_O)
+ $(CC) $< $(DSPL_O) -o $@ $(LFLAGS)
+
+$(OBJ_DIR)/%.o:$(SRC_DIR)/%.c
+ $(CC) $(CFLAGS) $< -o $@ $(LFLAGS)
+
+$(DSPL_O):$(DSPL_C)
+ $(CC) $(CFLAGS) $(DSPL_C) -o $(DSPL_O)
+
+bin/$(LIB_NAME):$(RELEASE_DIR)/$(LIB_NAME)
+ cp $(RELEASE_DIR)/$(LIB_NAME) bin/$(LIB_NAME)
+
+clean:
+ rm -f $(OBJ_DIR)/*.o
+ rm -f bin/*.exe
+
diff --git a/verification/bin/dat/.gitignore b/verification/bin/dat/.gitignore
new file mode 100644
index 0000000..b97d62b
--- /dev/null
+++ b/verification/bin/dat/.gitignore
@@ -0,0 +1,3 @@
+*.txt
+*.bin
+*.dat
\ No newline at end of file
diff --git a/verification/bin/octave/readbin.m b/verification/bin/octave/readbin.m
new file mode 100644
index 0000000..a37ddf0
--- /dev/null
+++ b/verification/bin/octave/readbin.m
@@ -0,0 +1,24 @@
+function [dat, n, m] = readbin(fn)
+
+ fid = fopen(fn);
+ if(~fid)
+ error('cannot to open file');
+ end
+
+ type = fread(fid, 1, 'int32');
+ n = fread(fid, 1, 'int32');
+ m = fread(fid, 1, 'int32');
+
+ if(type==0)
+ dat = fread(fid, [n*m, 1], 'double');
+ end
+
+ if(type==1)
+ y = fread(fid, [n*m*2, 1], 'double');
+ dat = y(1:2:end) + 1i * y(2:2:end);
+ end
+
+ dat = reshape(dat, n, m);
+
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/verification/bin/octave/writebin.m b/verification/bin/octave/writebin.m
new file mode 100644
index 0000000..c9c6e2b
--- /dev/null
+++ b/verification/bin/octave/writebin.m
@@ -0,0 +1,64 @@
+function res = writebin(x, type, fn)
+
+ if(type~=0 && type~=1)
+ res = 2;
+ return;
+ end
+
+ fid = fopen(fn, 'w');
+ if(~fid)
+ error('cannot to open file');
+ end
+
+ n = size(x, 1);
+ m = size(x, 2);
+
+ count = fwrite(fid, type, 'int32');
+ if(count ~= 1)
+ res = 1;
+ return;
+ end
+ count = fwrite(fid, n, 'int32');
+ if(count ~= 1)
+ fclose(fid);
+ res = 1;
+ return;
+ end
+ count = fwrite(fid, m, 'int32');
+ if(count ~= 1)
+ res = 1;
+ fclose(fid);
+ return;
+ end
+
+ flag = 0;
+ if(type==0)
+
+ count = fwrite(fid, x, 'double');
+ if(count ~= n*m)
+ res = 1;
+ fclose(fid);
+ return;
+ end
+ flag = 1;
+
+ else
+ y = reshape(x, n*m, 1);
+ z = zeros(2*n*m, 1);
+ z(1:2:end) = real(y);
+ z(2:2:end) = imag(y);
+ count = fwrite(fid, z, 'double');
+ if(count ~= 2*n*m)
+ res = 1;
+ fclose(fid);
+ return;
+ end
+ flag = 1;
+ end
+ if(flag == 0)
+ res = 3;
+ else
+ res = 0;
+ end
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/verification/bin/octave/writebin_readbin_verification.m b/verification/bin/octave/writebin_readbin_verification.m
new file mode 100644
index 0000000..7e74bdd
--- /dev/null
+++ b/verification/bin/octave/writebin_readbin_verification.m
@@ -0,0 +1,9 @@
+clear all; close all; clc;
+addpath('octave');
+
+x = readbin('dat/x.dat');
+if(isreal(x))
+ writebin(x, 0, 'dat/y.dat');
+else
+ writebin(x, 1, 'dat/y.dat');
+end
diff --git a/verification/bin/python/.gitignore b/verification/bin/python/.gitignore
new file mode 100644
index 0000000..e69de29
diff --git a/verification/obj/.gitignore b/verification/obj/.gitignore
new file mode 100644
index 0000000..444cfd9
--- /dev/null
+++ b/verification/obj/.gitignore
@@ -0,0 +1 @@
+.*o
\ No newline at end of file
diff --git a/verification/src/writebin_readbin_verification_complex.c b/verification/src/writebin_readbin_verification_complex.c
new file mode 100644
index 0000000..f7cd74e
--- /dev/null
+++ b/verification/src/writebin_readbin_verification_complex.c
@@ -0,0 +1,72 @@
+#include
+#include
+#include
+#include
+#include "dspl.h"
+
+int main(int argc, char* argv[])
+{
+ void* hdspl; /* DSPL handle */
+ complex_t *pxd = NULL;
+ complex_t *pyd = NULL;
+
+ int nx, ny, tx, err, verr;
+ double derr;
+
+ random_t rnd = {0}; /* random structure */
+
+ hdspl = dspl_load(); /* Load DSPL function */
+
+ printf("writebin and readbin for complex data:.........");
+
+ err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ err = randi(&nx, 1, 1, 1000000, &rnd);
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ pxd = (complex_t*) malloc(nx * sizeof(complex_t));
+
+ err = randn((double*)pxd, 2*nx, 1.0, 10.0, &rnd);
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ err = writebin((void*)pxd, nx, DAT_COMPLEX, "dat/x.dat");
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ err = system("octave octave/writebin_readbin_verification.m");
+
+ err = readbin("dat/y.dat", (void**)(&pyd), &ny, &tx);
+ if(err != RES_OK)
+ goto exit_error_code;
+ if(tx!=DAT_COMPLEX)
+ {
+ err = ERROR_DAT_TYPE;
+ goto exit_error_code;
+ }
+ verr = verif_cmplx(pxd, pyd, ny, 1E-12, &derr);
+ if(verr != DSPL_VERIF_SUCCESS)
+ goto exit_error_verif;
+
+ printf("ok (err = %12.4E)\n", derr);
+ goto exit_label;
+
+exit_error_code:
+ printf("FAILED (with code = 0x%8x)\n", err);
+ goto exit_label;
+exit_error_verif:
+ printf("FAILED (err = %12.4E)\n", derr);
+
+exit_label:
+ if(pxd)
+ free(pxd);
+ if(pyd)
+ free(pyd);
+ /* free dspl handle */
+ dspl_free(hdspl);
+ return 0;
+}
+
diff --git a/verification/src/writebin_readbin_verification_double.c b/verification/src/writebin_readbin_verification_double.c
new file mode 100644
index 0000000..7538209
--- /dev/null
+++ b/verification/src/writebin_readbin_verification_double.c
@@ -0,0 +1,73 @@
+#include
+#include
+#include
+#include
+#include "dspl.h"
+
+int main(int argc, char* argv[])
+{
+ void* hdspl; /* DSPL handle */
+ double *pxd = NULL;
+ double *pyd = NULL;
+
+ int nx, ny, tx, err, verr;
+ double derr;
+
+ random_t rnd = {0}; /* random structure */
+
+ hdspl = dspl_load(); /* Load DSPL function */
+
+ printf("writebin and readbin for double data:..........");
+
+ err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ err = randi(&nx, 1, 1, 1000000, &rnd);
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ pxd = (double*) malloc(nx * sizeof(double));
+
+ err = randn(pxd, nx, 1.0, 10.0, &rnd);
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ err = writebin((void*)pxd, nx, DAT_DOUBLE, "dat/x.dat");
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ err = system("octave octave/writebin_readbin_verification.m");
+
+ err = readbin("dat/y.dat", (void**)(&pyd), &ny, &tx);
+ if(err != RES_OK)
+ goto exit_error_code;
+
+ if(tx!=DAT_DOUBLE)
+ {
+ err = ERROR_DAT_TYPE;
+ goto exit_error_code;
+ }
+ verr = verif(pxd, pyd, ny, 1E-12, &derr);
+ if(verr != DSPL_VERIF_SUCCESS)
+ goto exit_error_verif;
+
+ printf("ok (err = %12.4E)", derr);
+ goto exit_label;
+
+exit_error_code:
+ printf("FAILED (with code = 0x%8x)", err);
+ goto exit_label;
+exit_error_verif:
+ printf("FAILED (err = %12.4E)", derr);
+
+exit_label:
+ if(pxd)
+ free(pxd);
+ if(pyd)
+ free(pyd);
+ /* free dspl handle */
+ dspl_free(hdspl);
+ return 0;
+}
+
diff --git a/verification/verif.sh b/verification/verif.sh
new file mode 100644
index 0000000..89d3243
--- /dev/null
+++ b/verification/verif.sh
@@ -0,0 +1,7 @@
+mingw32-make
+
+cd bin
+for file in *.exe
+do
+ "./$file"
+done