kopia lustrzana https://github.com/Dsplib/libdspl-2.0
1) added xcorr fucntions for real and complex data
2) added performance test direcories 3) added verification directories and files for readbin abd writebin verification with octave Changes to be committed: modified: Makefile modified: dspl/src/inout.c modified: dspl/src/statistic.c new file: dspl/src/xcorr.c new file: examples/src/xcorr_cmplx_test.c new file: examples/src/xcorr_test.c modified: include/dspl.c modified: include/dspl.h new file: performance/Makefile new file: performance/bin/.gitignore new file: performance/obj/.gitignore new file: performance/src/xcorr_cmplx_performance.c new file: verification/Makefile new file: verification/bin/dat/.gitignore new file: verification/bin/octave/readbin.m new file: verification/bin/octave/writebin.m new file: verification/bin/octave/writebin_readbin_verification.m new file: verification/bin/python/.gitignore new file: verification/obj/.gitignore new file: verification/src/writebin_readbin_verification_complex.c new file: verification/src/writebin_readbin_verification_double.c new file: verification/verif.shpull/6/merge
rodzic
52a33913b3
commit
23c2496b2a
7
Makefile
7
Makefile
|
@ -6,8 +6,11 @@ all:
|
||||||
$(MAKE) -C lapack
|
$(MAKE) -C lapack
|
||||||
$(MAKE) -C dspl
|
$(MAKE) -C dspl
|
||||||
$(MAKE) -C examples
|
$(MAKE) -C examples
|
||||||
|
$(MAKE) -C performance
|
||||||
|
$(MAKE) -C verification
|
||||||
|
|
||||||
clean:
|
clean:
|
||||||
$(MAKE) -C dspl clean
|
$(MAKE) -C dspl clean
|
||||||
$(MAKE) -C examples clean
|
$(MAKE) -C examples clean
|
||||||
|
$(MAKE) -C performance clean
|
||||||
|
$(MAKE) -C verification clean
|
|
@ -44,9 +44,6 @@ void DSPL_API dspl_info()
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#ifdef DOXYGEN_ENGLISH
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
@ -90,7 +87,6 @@ int DSPL_API readbin(char* fn, void** x, int* k, int* dtype)
|
||||||
|
|
||||||
if(k)
|
if(k)
|
||||||
*k = n*m;
|
*k = n*m;
|
||||||
|
|
||||||
switch(t)
|
switch(t)
|
||||||
{
|
{
|
||||||
case DAT_DOUBLE:
|
case DAT_DOUBLE:
|
||||||
|
@ -111,6 +107,7 @@ int DSPL_API readbin(char* fn, void** x, int* k, int* dtype)
|
||||||
err = ERROR_FREAD_SIZE;
|
err = ERROR_FREAD_SIZE;
|
||||||
goto exit_label;
|
goto exit_label;
|
||||||
}
|
}
|
||||||
|
break;
|
||||||
default:
|
default:
|
||||||
err = ERROR_DAT_TYPE;
|
err = ERROR_DAT_TYPE;
|
||||||
goto exit_label;
|
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)
|
if(fwrite(&n, sizeof(int), 1, pFile) != 1)
|
||||||
{
|
{
|
||||||
res = ERROR_FWRITE_SIZE;
|
res = ERROR_FWRITE_SIZE;
|
||||||
|
|
|
@ -257,4 +257,8 @@ int DSPL_API minmax(double* x, int n, double* xmin, double* xmax)
|
||||||
*xmax = max;
|
*xmax = max;
|
||||||
|
|
||||||
return RES_OK;
|
return RES_OK;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#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;
|
||||||
|
}
|
|
@ -0,0 +1,78 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#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;
|
||||||
|
}
|
||||||
|
|
|
@ -0,0 +1,63 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#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;
|
||||||
|
}
|
||||||
|
|
|
@ -185,6 +185,9 @@ p_writetxt_cmplx_im writetxt_cmplx_im ;
|
||||||
p_writetxt_cmplx_re writetxt_cmplx_re ;
|
p_writetxt_cmplx_re writetxt_cmplx_re ;
|
||||||
p_writetxt_int writetxt_int ;
|
p_writetxt_int writetxt_int ;
|
||||||
|
|
||||||
|
p_xcorr xcorr ;
|
||||||
|
p_xcorr_cmplx xcorr_cmplx ;
|
||||||
|
|
||||||
|
|
||||||
#ifdef WIN_OS
|
#ifdef WIN_OS
|
||||||
#define LOAD_FUNC(fn) \
|
#define LOAD_FUNC(fn) \
|
||||||
|
@ -488,6 +491,10 @@ void* dspl_load()
|
||||||
LOAD_FUNC(writetxt_cmplx_re);
|
LOAD_FUNC(writetxt_cmplx_re);
|
||||||
LOAD_FUNC(writetxt_int);
|
LOAD_FUNC(writetxt_int);
|
||||||
|
|
||||||
|
LOAD_FUNC(xcorr);
|
||||||
|
LOAD_FUNC(xcorr_cmplx);
|
||||||
|
|
||||||
|
|
||||||
#ifdef WIN_OS
|
#ifdef WIN_OS
|
||||||
return (void*)handle;
|
return (void*)handle;
|
||||||
exit_label:
|
exit_label:
|
||||||
|
|
|
@ -561,6 +561,7 @@ Variable `z = 1-2j`, here `j` - imaginary unit, but variable `y = 5`.
|
||||||
#define ERROR_WIN_SYM 0x23091925
|
#define ERROR_WIN_SYM 0x23091925
|
||||||
#define ERROR_WIN_TYPE 0x23092025
|
#define ERROR_WIN_TYPE 0x23092025
|
||||||
/* X 0x24xxxxxx*/
|
/* X 0x24xxxxxx*/
|
||||||
|
#define ERROR_XCORR_FLAG 0x24031518
|
||||||
/* Y 0x25xxxxxx*/
|
/* Y 0x25xxxxxx*/
|
||||||
/* Z 0x26xxxxxx*/
|
/* 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_FILTER_ELLIP 0x00000800
|
||||||
|
|
||||||
|
|
||||||
|
#define DSPL_XCORR_NOSCALE 0x00000000
|
||||||
|
#define DSPL_XCORR_BIASED 0x00000001
|
||||||
|
#define DSPL_XCORR_UNBIASED 0x00000002
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#define ELLIP_ITER 16
|
#define ELLIP_ITER 16
|
||||||
#define ELLIP_MAX_ORD 24
|
#define ELLIP_MAX_ORD 24
|
||||||
|
@ -1400,6 +1406,24 @@ DECLARE_FUNC(int, writetxt_int, int*
|
||||||
COMMA int
|
COMMA int
|
||||||
COMMA char*);
|
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
|
#ifdef __cplusplus
|
||||||
|
|
|
@ -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
|
||||||
|
|
|
@ -0,0 +1,2 @@
|
||||||
|
*.exe
|
||||||
|
*.dll
|
|
@ -0,0 +1 @@
|
||||||
|
.*o
|
|
@ -0,0 +1,57 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <time.h>
|
||||||
|
#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;
|
||||||
|
}
|
||||||
|
|
|
@ -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
|
||||||
|
|
|
@ -0,0 +1,3 @@
|
||||||
|
*.txt
|
||||||
|
*.bin
|
||||||
|
*.dat
|
|
@ -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
|
|
@ -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
|
|
@ -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
|
|
@ -0,0 +1 @@
|
||||||
|
.*o
|
|
@ -0,0 +1,72 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <time.h>
|
||||||
|
#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;
|
||||||
|
}
|
||||||
|
|
|
@ -0,0 +1,73 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <time.h>
|
||||||
|
#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;
|
||||||
|
}
|
||||||
|
|
|
@ -0,0 +1,7 @@
|
||||||
|
mingw32-make
|
||||||
|
|
||||||
|
cd bin
|
||||||
|
for file in *.exe
|
||||||
|
do
|
||||||
|
"./$file"
|
||||||
|
done
|
Ładowanie…
Reference in New Issue