kopia lustrzana https://github.com/Dsplib/libdspl-2.0
doxyfiles has changed;
added performance estimation for FFT added verificator for FFT and statistic functions (mean and std) Changes to be committed: modified: dspl/dox/doxyfile_en modified: dspl/dox/doxyfile_ru new file: dspl/dox/doxygen_equation.css modified: dspl/dox/header_ru.html modified: dspl/dox/ru/mainpage.dox modified: dspl/src/array.c modified: dspl/src/fft.c modified: dspl/src/fft_subkernel.c modified: dspl/src/filter_an.c modified: dspl/src/inout.c modified: dspl/src/statistic.c new file: dspl/src/verification.c modified: dspl/src/xcorr.c modified: include/dspl.c modified: include/dspl.h new file: performance/bin/dat/.gitignore new file: performance/bin/octave/fft_cmplx_performance.m new file: performance/bin/python/fft_cmplx_performance.py new file: performance/src/fft_cmplx_performance.c modified: verification/Makefile new file: verification/bin/octave/fft_verification.m new file: verification/bin/octave/statistic_verification.m modified: verification/bin/octave/writebin_readbin_verification.m deleted: verification/src/dspl_verif.h new file: verification/src/fft_verification.c new file: verification/src/statistic_verification.c new file: verification/src/writebin_readbin_verification.c deleted: verification/src/writebin_readbin_verification_complex.c deleted: verification/src/writebin_readbin_verification_double.cpull/6/merge
rodzic
226b1ff7d5
commit
d2e9224f37
|
@ -1551,7 +1551,7 @@ GENERATE_TREEVIEW = NO
|
||||||
# Minimum value: 0, maximum value: 20, default value: 4.
|
# Minimum value: 0, maximum value: 20, default value: 4.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
ENUM_VALUES_PER_LINE = 8
|
ENUM_VALUES_PER_LINE = 4
|
||||||
|
|
||||||
# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be used
|
# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be used
|
||||||
# to set the initial width (in pixels) of the frame in which the tree is shown.
|
# to set the initial width (in pixels) of the frame in which the tree is shown.
|
||||||
|
@ -1576,7 +1576,7 @@ EXT_LINKS_IN_WINDOW = NO
|
||||||
# The default value is: png.
|
# The default value is: png.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
HTML_FORMULA_FORMAT = png
|
HTML_FORMULA_FORMAT = svg
|
||||||
|
|
||||||
# Use this tag to change the font size of LaTeX formulas included as images in
|
# Use this tag to change the font size of LaTeX formulas included as images in
|
||||||
# the HTML documentation. When you change the font size after a successful
|
# the HTML documentation. When you change the font size after a successful
|
||||||
|
@ -1585,7 +1585,7 @@ HTML_FORMULA_FORMAT = png
|
||||||
# Minimum value: 8, maximum value: 50, default value: 10.
|
# Minimum value: 8, maximum value: 50, default value: 10.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
FORMULA_FONTSIZE = 14
|
FORMULA_FONTSIZE = 12
|
||||||
|
|
||||||
# Use the FORMULA_TRANSPARENT tag to determine whether or not the images
|
# Use the FORMULA_TRANSPARENT tag to determine whether or not the images
|
||||||
# generated for formulas are transparent PNGs. Transparent PNGs are not
|
# generated for formulas are transparent PNGs. Transparent PNGs are not
|
||||||
|
@ -1596,7 +1596,7 @@ FORMULA_FONTSIZE = 14
|
||||||
# The default value is: YES.
|
# The default value is: YES.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
FORMULA_TRANSPARENT = YES
|
FORMULA_TRANSPARENT = NO
|
||||||
|
|
||||||
# The FORMULA_MACROFILE can contain LaTeX \newcommand and \renewcommand commands
|
# The FORMULA_MACROFILE can contain LaTeX \newcommand and \renewcommand commands
|
||||||
# to create new LaTeX commands to be used in formulas as building blocks. See
|
# to create new LaTeX commands to be used in formulas as building blocks. See
|
||||||
|
@ -1613,7 +1613,7 @@ FORMULA_MACROFILE =
|
||||||
# The default value is: NO.
|
# The default value is: NO.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
USE_MATHJAX = YES
|
USE_MATHJAX = NO
|
||||||
|
|
||||||
# When MathJax is enabled you can set the default output format to be used for
|
# When MathJax is enabled you can set the default output format to be used for
|
||||||
# the MathJax output. See the MathJax site (see:
|
# the MathJax output. See the MathJax site (see:
|
||||||
|
@ -1623,7 +1623,7 @@ USE_MATHJAX = YES
|
||||||
# The default value is: HTML-CSS.
|
# The default value is: HTML-CSS.
|
||||||
# This tag requires that the tag USE_MATHJAX is set to YES.
|
# This tag requires that the tag USE_MATHJAX is set to YES.
|
||||||
|
|
||||||
MATHJAX_FORMAT = HTML-CSS
|
MATHJAX_FORMAT = SVG
|
||||||
|
|
||||||
# When MathJax is enabled you need to specify the location relative to the HTML
|
# When MathJax is enabled you need to specify the location relative to the HTML
|
||||||
# output directory using the MATHJAX_RELPATH option. The destination directory
|
# output directory using the MATHJAX_RELPATH option. The destination directory
|
||||||
|
@ -1636,7 +1636,7 @@ MATHJAX_FORMAT = HTML-CSS
|
||||||
# The default value is: https://cdn.jsdelivr.net/npm/mathjax@2.
|
# The default value is: https://cdn.jsdelivr.net/npm/mathjax@2.
|
||||||
# This tag requires that the tag USE_MATHJAX is set to YES.
|
# This tag requires that the tag USE_MATHJAX is set to YES.
|
||||||
|
|
||||||
MATHJAX_RELPATH = http://cdn.mathjax.org/mathjax/latest
|
MATHJAX_RELPATH = http://dsplib.org/mathjax/latest
|
||||||
|
|
||||||
# The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax
|
# The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax
|
||||||
# extension names that should be enabled during MathJax rendering. For example
|
# extension names that should be enabled during MathJax rendering. For example
|
||||||
|
|
|
@ -1551,7 +1551,7 @@ GENERATE_TREEVIEW = NO
|
||||||
# Minimum value: 0, maximum value: 20, default value: 4.
|
# Minimum value: 0, maximum value: 20, default value: 4.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
ENUM_VALUES_PER_LINE = 8
|
ENUM_VALUES_PER_LINE = 4
|
||||||
|
|
||||||
# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be used
|
# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be used
|
||||||
# to set the initial width (in pixels) of the frame in which the tree is shown.
|
# to set the initial width (in pixels) of the frame in which the tree is shown.
|
||||||
|
@ -1576,7 +1576,7 @@ EXT_LINKS_IN_WINDOW = NO
|
||||||
# The default value is: png.
|
# The default value is: png.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
HTML_FORMULA_FORMAT = png
|
HTML_FORMULA_FORMAT = svg
|
||||||
|
|
||||||
# Use this tag to change the font size of LaTeX formulas included as images in
|
# Use this tag to change the font size of LaTeX formulas included as images in
|
||||||
# the HTML documentation. When you change the font size after a successful
|
# the HTML documentation. When you change the font size after a successful
|
||||||
|
@ -1585,7 +1585,7 @@ HTML_FORMULA_FORMAT = png
|
||||||
# Minimum value: 8, maximum value: 50, default value: 10.
|
# Minimum value: 8, maximum value: 50, default value: 10.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
FORMULA_FONTSIZE = 14
|
FORMULA_FONTSIZE = 12
|
||||||
|
|
||||||
# Use the FORMULA_TRANSPARENT tag to determine whether or not the images
|
# Use the FORMULA_TRANSPARENT tag to determine whether or not the images
|
||||||
# generated for formulas are transparent PNGs. Transparent PNGs are not
|
# generated for formulas are transparent PNGs. Transparent PNGs are not
|
||||||
|
@ -1613,7 +1613,7 @@ FORMULA_MACROFILE =
|
||||||
# The default value is: NO.
|
# The default value is: NO.
|
||||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||||
|
|
||||||
USE_MATHJAX = YES
|
USE_MATHJAX = NO
|
||||||
|
|
||||||
# When MathJax is enabled you can set the default output format to be used for
|
# When MathJax is enabled you can set the default output format to be used for
|
||||||
# the MathJax output. See the MathJax site (see:
|
# the MathJax output. See the MathJax site (see:
|
||||||
|
|
|
@ -0,0 +1,5 @@
|
||||||
|
.formulaInl {
|
||||||
|
-ms-transform: scale(1.5, 1.5); /* IE 9 */
|
||||||
|
-webkit-transform: scale(1.5, 1.5); /* Chrome, Safari, Opera */
|
||||||
|
transform: scale(1.5, 1.5);
|
||||||
|
}
|
|
@ -17,6 +17,10 @@
|
||||||
<link href="$relpath^tabs.css" rel="stylesheet" type="text/css"/>
|
<link href="$relpath^tabs.css" rel="stylesheet" type="text/css"/>
|
||||||
<script type="text/javascript" src="$relpath^jquery.js"></script>
|
<script type="text/javascript" src="$relpath^jquery.js"></script>
|
||||||
<script type="text/javascript" src="$relpath^dynsections.js"></script>
|
<script type="text/javascript" src="$relpath^dynsections.js"></script>
|
||||||
|
|
||||||
|
<!-- Place this tag in your head or just before your close body tag. -->
|
||||||
|
<script async defer src="https://buttons.github.io/buttons.js"></script>
|
||||||
|
|
||||||
$treeview
|
$treeview
|
||||||
$search
|
$search
|
||||||
$mathjax
|
$mathjax
|
||||||
|
|
|
@ -13,11 +13,9 @@ DSPL-2.0 --- свободная библиотека алгоритмов циф
|
||||||
условии динамической линковки.
|
условии динамической линковки.
|
||||||
|
|
||||||
|
|
||||||
Исходные коды библиотеки доступны на
|
<!-- Place this tag where you want the button to render. -->
|
||||||
<a href = "https://github.com/Dsplib/libdspl-2.0">GitHub</a>. \n
|
<a class="github-button" href="https://github.com/Dsplib/libdspl-2.0" data-color-scheme="no-preference: light; light: light; dark: light;" data-size="large" data-show-count="true" aria-label="Follow @Dsplib on GitHub">Исходные коды библиотеки libdspl-2.0 на GitHub</a>
|
||||||
|
|
||||||
Вы также можете внести свой вклад в развитие данной библиотеки.
|
|
||||||
\b Присоединяйтесь!
|
|
||||||
|
|
||||||
|
|
||||||
\subsection sec_doc_start Компиляция и запуск программ DSPL-2.0
|
\subsection sec_doc_start Компиляция и запуск программ DSPL-2.0
|
||||||
|
|
291
dspl/src/array.c
291
dspl/src/array.c
|
@ -1310,294 +1310,3 @@ int DSPL_API ones(double* x, int n)
|
||||||
x[i] = 1.0;
|
x[i] = 1.0;
|
||||||
return RES_OK;
|
return RES_OK;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
#ifdef DOXYGEN_ENGLISH
|
|
||||||
/*! ****************************************************************************
|
|
||||||
\ingroup ARRAY_GROUP
|
|
||||||
\fn int verif(double* x, double* y, size_t n, double eps, double* err)
|
|
||||||
\brief Real arrays verification
|
|
||||||
|
|
||||||
Function calculates a maximum relative error between two real arrays `x`
|
|
||||||
and `y` (both length equals `n`):
|
|
||||||
|
|
||||||
\f[
|
|
||||||
e = \max \left( \frac{|x(k) - y(k)| }{ |x(k)|} \right),
|
|
||||||
\quad if \quad |x(k)| > 0,
|
|
||||||
\f]
|
|
||||||
or
|
|
||||||
\f[
|
|
||||||
e = \max(|x(k) - y(k)| ), ~\qquad if \quad~|x(k)| = 0,
|
|
||||||
\f]
|
|
||||||
|
|
||||||
This function can be used for algorithms verification if vector `x` is user
|
|
||||||
algorithm result and vector `y` -- reference vector.
|
|
||||||
|
|
||||||
\param[in] x
|
|
||||||
Pointer to the first vector `x`. \n
|
|
||||||
Vector size is `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] y
|
|
||||||
Pointer to the second vector `y`. \n
|
|
||||||
Vector size is `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] n
|
|
||||||
Size of vectors `x` and `y`. \n \n
|
|
||||||
|
|
||||||
\param[in] eps
|
|
||||||
Relative error threshold. \n
|
|
||||||
If error less than `eps`, then function returns
|
|
||||||
`DSPL_VERIF_SUCCESS`, else `DSPL_VERIF_FAILED`. \n \n
|
|
||||||
|
|
||||||
\param[in, out] err
|
|
||||||
Pointer to the variable which keep
|
|
||||||
maximum relative error. \n
|
|
||||||
Pointer can be `NULL`, maximum error will not be returned
|
|
||||||
in this case. \n \n
|
|
||||||
|
|
||||||
\return
|
|
||||||
`DSPL_VERIF_SUCCESS` if maximum relative error less than `eps`. \n
|
|
||||||
Otherwise `DSPL_VERIF_FAILED`.
|
|
||||||
|
|
||||||
\author Sergey Bakhurin www.dsplib.org
|
|
||||||
***************************************************************************** */
|
|
||||||
#endif
|
|
||||||
#ifdef DOXYGEN_RUSSIAN
|
|
||||||
/*! ****************************************************************************
|
|
||||||
\ingroup ARRAY_GROUP
|
|
||||||
\fn int verif(double* x, double* y, size_t n, double eps, double* err)
|
|
||||||
\brief Верификация вещественных массивов
|
|
||||||
|
|
||||||
Функция производит расчет максимальной относительной ошибки между вещественными
|
|
||||||
векторами `x` и `y` равной длины `n`:
|
|
||||||
|
|
||||||
\f[
|
|
||||||
e = \max \left( \frac{|x(k) - y(k)| }{ |x(k)|} \right), \quad \quad |x(k)| > 0,
|
|
||||||
\f]
|
|
||||||
или
|
|
||||||
\f[
|
|
||||||
e = \max(|x(k) - y(k)| ), \qquad \quad~|x(k)| = 0,
|
|
||||||
\f]
|
|
||||||
и возвращает `DSPL_VERIF_SUCCESS` если
|
|
||||||
разница \f$ e\f$ меньше `eps`.
|
|
||||||
В противном случае возвращает `DSPL_VERIF_FAILED`. \n
|
|
||||||
Данная функция используется для верификации работы алгоритмов если вектор `x`
|
|
||||||
результат работы алгоритма пользователя, а `y` -- результат работы этого же
|
|
||||||
алгоритма сторонней функцией.
|
|
||||||
|
|
||||||
\param[in] x
|
|
||||||
Указатель на первый вектор `x`. \n
|
|
||||||
Размер вектора `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] y
|
|
||||||
Указатель на второй вектор `y`. \n
|
|
||||||
Размер вектора `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] n
|
|
||||||
Размер векторов `x` и `y`. \n \n
|
|
||||||
|
|
||||||
\param[in] eps
|
|
||||||
Допустимая относительная ошибка. \n
|
|
||||||
Если максимальная относительная ошибка меньше `eps`, то функция возвращает
|
|
||||||
`DSPL_VERIF_SUCCESS`, в противном случае `DSPL_VERIF_FAILED`. \n \n
|
|
||||||
|
|
||||||
\param[in, out] err
|
|
||||||
Указатель на переменную максимальной относительной ошибки. \n
|
|
||||||
По данному адресу будет записано значение максимальной относительной ошибки. \n
|
|
||||||
Указатель может быть `NULL`, значение ошибки в этом случае
|
|
||||||
не возвращается. \n \n
|
|
||||||
|
|
||||||
\return
|
|
||||||
`DSPL_VERIF_SUCCESS` если относительная ошибка меньше `eps`. \n
|
|
||||||
В противном случае `DSPL_VERIF_FAILED`.
|
|
||||||
|
|
||||||
\author Бахурин Сергей www.dsplib.org
|
|
||||||
***************************************************************************** */
|
|
||||||
#endif
|
|
||||||
int DSPL_API verif(double* x, double* y, size_t n, double eps, double* err)
|
|
||||||
{
|
|
||||||
double d, maxd;
|
|
||||||
size_t k;
|
|
||||||
int res;
|
|
||||||
if(!x || !y)
|
|
||||||
return ERROR_PTR;
|
|
||||||
if(n < 1)
|
|
||||||
return ERROR_SIZE;
|
|
||||||
if(eps <= 0.0 )
|
|
||||||
return ERROR_NEGATIVE;
|
|
||||||
|
|
||||||
maxd = -100.0;
|
|
||||||
|
|
||||||
for(k = 0; k < n; k++)
|
|
||||||
{
|
|
||||||
d = fabs(x[k] - y[k]);
|
|
||||||
if(fabs(x[k]) > 0.0)
|
|
||||||
{
|
|
||||||
d = d / fabs(x[k]);
|
|
||||||
if(d > maxd)
|
|
||||||
maxd = d;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if(err)
|
|
||||||
*err = maxd;
|
|
||||||
|
|
||||||
if(maxd > eps)
|
|
||||||
res = DSPL_VERIF_FAILED;
|
|
||||||
else
|
|
||||||
res = DSPL_VERIF_SUCCESS;
|
|
||||||
|
|
||||||
return res;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#ifdef DOXYGEN_ENGLISH
|
|
||||||
/*! ****************************************************************************
|
|
||||||
\ingroup ARRAY_GROUP
|
|
||||||
\fn int verif_cmplx(complex_t* x, complex_t* y, size_t n,
|
|
||||||
double eps, double* err)
|
|
||||||
\brief
|
|
||||||
Complex arrays verification
|
|
||||||
|
|
||||||
Function calculates a maximum relative error between two complex arrays `x`
|
|
||||||
and `y` (both length equals `n`):
|
|
||||||
|
|
||||||
\f[
|
|
||||||
e = \max \left( \frac{|x(k) - y(k)| }{ |x(k)|} \right),
|
|
||||||
\quad if \quad |x(k)| > 0,
|
|
||||||
\f]
|
|
||||||
or
|
|
||||||
\f[
|
|
||||||
e = \max(|x(k) - y(k)| ), ~\qquad if \quad~|x(k)| = 0,
|
|
||||||
\f]
|
|
||||||
Return `DSPL_VERIF_SUCCESS` if maximum relative error \f$ e\f$ less than `eps`.
|
|
||||||
Else returns `DSPL_VERIF_FAILED`. \n
|
|
||||||
|
|
||||||
This function can be used for algorithms verification if vector `x` is user
|
|
||||||
algorithm result and vector `y` -- reference vector.
|
|
||||||
|
|
||||||
\param[in] x
|
|
||||||
Pointer to the first vector `x`. \n
|
|
||||||
Vector size is `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] y
|
|
||||||
Pointer to the second vector `y`. \n
|
|
||||||
Vector size is `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] n
|
|
||||||
Size of vectors `x` and `y`. \n \n
|
|
||||||
|
|
||||||
\param[in] eps
|
|
||||||
Relative error threshold. \n
|
|
||||||
If error less than `eps`, then function returns
|
|
||||||
`DSPL_VERIF_SUCCESS`, else `DSPL_VERIF_FAILED`. \n \n
|
|
||||||
|
|
||||||
\param[in, out] err
|
|
||||||
Pointer to the variable which keep
|
|
||||||
maximum relative error. \n
|
|
||||||
Pointer can be `NULL`, maximum error will not be returned
|
|
||||||
in this case. \n \n
|
|
||||||
|
|
||||||
\return
|
|
||||||
`DSPL_VERIF_SUCCESS` if maximum relative error less than `eps`. \n
|
|
||||||
Otherwise `DSPL_VERIF_FAILED`.
|
|
||||||
|
|
||||||
\author
|
|
||||||
Sergey Bakhurin
|
|
||||||
www.dsplib.org
|
|
||||||
***************************************************************************** */
|
|
||||||
#endif
|
|
||||||
#ifdef DOXYGEN_RUSSIAN
|
|
||||||
/*! ****************************************************************************
|
|
||||||
\ingroup ARRAY_GROUP
|
|
||||||
\fn int verif_cmplx(complex_t* x, complex_t* y, size_t n,
|
|
||||||
double eps, double* err)
|
|
||||||
\brief Верификация комплексных массивов
|
|
||||||
|
|
||||||
Функция производит расчет максимальной относительной ошибки между комплексными
|
|
||||||
векторами `x` и `y` равной длины `n`:
|
|
||||||
|
|
||||||
\f[
|
|
||||||
e = \max \left( \frac{|x(k) - y(k)|}{|x(k)|} \right), \quad \quad |x(k)| > 0,
|
|
||||||
\f]
|
|
||||||
или
|
|
||||||
\f[
|
|
||||||
e = \max(|x(k) - y(k)| ), ~\qquad \quad~|x(k)| = 0,
|
|
||||||
\f]
|
|
||||||
и возвращает `DSPL_VERIF_SUCCESS` если
|
|
||||||
разница \f$ e\f$ меньше `eps`.
|
|
||||||
В противном случае возвращает `DSPL_VERIF_FAILED`. \n
|
|
||||||
Данная функция используется для верификации работы алгоритмов если вектор `x`
|
|
||||||
результат работы алгоритма пользователя, а `y` -- результат работы этого же
|
|
||||||
алгоритма сторонней функцией.
|
|
||||||
|
|
||||||
\param[in] x
|
|
||||||
Указатель на первый вектор `x`. \n
|
|
||||||
Размер вектора `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] y
|
|
||||||
Указатель на второй вектор `y`. \n
|
|
||||||
Размер вектора `[n x 1]`. \n \n
|
|
||||||
|
|
||||||
\param[in] n
|
|
||||||
Размер векторов `x` и `y`. \n \n
|
|
||||||
|
|
||||||
\param[in] eps
|
|
||||||
Допустимая относительная ошибка. \n
|
|
||||||
Если максимальная относительная ошибка меньше `eps`, то функция возвращает
|
|
||||||
`DSPL_VERIF_SUCCESS`, в противном случае `DSPL_VERIF_FAILED`. \n \n
|
|
||||||
|
|
||||||
\param[in, out] err
|
|
||||||
Указатель на переменную максимальной относительной ошибки. \n
|
|
||||||
По данному адресу будет записано значение максимальной относительной ошибки. \n
|
|
||||||
Указатель может быть `NULL`, значение ошибки в этом
|
|
||||||
случае не возвращается. \n \n
|
|
||||||
|
|
||||||
\return
|
|
||||||
`DSPL_VERIF_SUCCESS` если функция выполнена успешно. \n
|
|
||||||
В противном случае `DSPL_VERIF_FAILED`.
|
|
||||||
|
|
||||||
\author
|
|
||||||
Бахурин Сергей www.dsplib.org
|
|
||||||
***************************************************************************** */
|
|
||||||
#endif
|
|
||||||
int DSPL_API verif_cmplx(complex_t* x, complex_t* y, size_t n,
|
|
||||||
double eps, double* err)
|
|
||||||
{
|
|
||||||
|
|
||||||
complex_t d;
|
|
||||||
double mx, md, maxd;
|
|
||||||
size_t k;
|
|
||||||
int res;
|
|
||||||
if(!x || !y)
|
|
||||||
return ERROR_PTR;
|
|
||||||
if(n < 1)
|
|
||||||
return ERROR_SIZE;
|
|
||||||
if(eps <= 0.0 )
|
|
||||||
return ERROR_NEGATIVE;
|
|
||||||
|
|
||||||
maxd = -100.0;
|
|
||||||
|
|
||||||
for(k = 0; k < n; k++)
|
|
||||||
{
|
|
||||||
RE(d) = RE(x[k]) - RE(y[k]);
|
|
||||||
IM(d) = IM(x[k]) - IM(y[k]);
|
|
||||||
md = ABS(d);
|
|
||||||
mx = ABS(x[k]);
|
|
||||||
if(mx > 0.0)
|
|
||||||
{
|
|
||||||
md = md / mx;
|
|
||||||
if(md > maxd)
|
|
||||||
maxd = md;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if(err)
|
|
||||||
*err = maxd;
|
|
||||||
|
|
||||||
if(maxd > eps)
|
|
||||||
res = DSPL_VERIF_FAILED;
|
|
||||||
else
|
|
||||||
res = DSPL_VERIF_SUCCESS;
|
|
||||||
|
|
||||||
return res;
|
|
||||||
}
|
|
||||||
|
|
|
@ -590,7 +590,7 @@ label_size:
|
||||||
memcpy(t1, t0, n*sizeof(complex_t));
|
memcpy(t1, t0, n*sizeof(complex_t));
|
||||||
matrix_transpose_cmplx(t1, n2, n1, t0);
|
matrix_transpose_cmplx(t1, n2, n1, t0);
|
||||||
}
|
}
|
||||||
|
|
||||||
if(n1 == 16)
|
if(n1 == 16)
|
||||||
for(k = 0; k < n2; k++)
|
for(k = 0; k < n2; k++)
|
||||||
dft16(t0+16*k, t1+16*k);
|
dft16(t0+16*k, t1+16*k);
|
||||||
|
@ -785,13 +785,13 @@ int DSPL_API fft_create(fft_t* pfft, int n)
|
||||||
while(s > 1)
|
while(s > 1)
|
||||||
{
|
{
|
||||||
n2 = 1;
|
n2 = 1;
|
||||||
if(s%16== 0) { n2 = 16; goto label_size; }
|
if(s%16== 0) { n2 = 16; goto label_size; }
|
||||||
if(s%7 == 0) { n2 = 7; goto label_size; }
|
if(s%7 == 0) { n2 = 7; goto label_size; }
|
||||||
if(s%8 == 0) { n2 = 8; goto label_size; }
|
if(s%8 == 0) { n2 = 8; goto label_size; }
|
||||||
if(s%5 == 0) { n2 = 5; goto label_size; }
|
if(s%5 == 0) { n2 = 5; goto label_size; }
|
||||||
if(s%4 == 0) { n2 = 4; goto label_size; }
|
if(s%4 == 0) { n2 = 4; goto label_size; }
|
||||||
if(s%3 == 0) { n2 = 3; goto label_size; }
|
if(s%3 == 0) { n2 = 3; goto label_size; }
|
||||||
if(s%2 == 0) { n2 = 2; goto label_size; }
|
if(s%2 == 0) { n2 = 2; goto label_size; }
|
||||||
|
|
||||||
|
|
||||||
label_size:
|
label_size:
|
||||||
|
|
|
@ -26,6 +26,10 @@
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/*******************************************************************************
|
/*******************************************************************************
|
||||||
2 points DFT
|
2 points DFT
|
||||||
*******************************************************************************/
|
*******************************************************************************/
|
||||||
|
@ -376,7 +380,7 @@ void dft8 (complex_t *x, complex_t* y)
|
||||||
/*******************************************************************************
|
/*******************************************************************************
|
||||||
16 points DFT (Winograd algorithm)
|
16 points DFT (Winograd algorithm)
|
||||||
*******************************************************************************/
|
*******************************************************************************/
|
||||||
void dft16 (complex_t *x, complex_t* y)
|
void dft16(complex_t *x, complex_t* y)
|
||||||
{
|
{
|
||||||
complex_t t0[16];
|
complex_t t0[16];
|
||||||
complex_t t1[16];
|
complex_t t1[16];
|
||||||
|
@ -451,8 +455,6 @@ void dft16 (complex_t *x, complex_t* y)
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/*******************************************************************************
|
/*******************************************************************************
|
||||||
4 x 2 matrix transpose
|
4 x 2 matrix transpose
|
||||||
*******************************************************************************/
|
*******************************************************************************/
|
||||||
|
|
|
@ -1064,7 +1064,7 @@ Phase delay calculation for digital or analog filter corresponds to
|
||||||
|
|
||||||
Group delay is describes as:
|
Group delay is describes as:
|
||||||
\f[
|
\f[
|
||||||
\tau_\{\varphi}(\omega) = - \frac{\Phi(\omega)}{\omega},
|
\tau_{\varphi}(\omega) = - \frac{\Phi(\omega)}{\omega},
|
||||||
\f]
|
\f]
|
||||||
here \f$\Phi(\omega)\f$ -- filter phase response, \f$\omega\f$ is angular
|
here \f$\Phi(\omega)\f$ -- filter phase response, \f$\omega\f$ is angular
|
||||||
frequency for analog filter, or normalized frequency for digital filter.
|
frequency for analog filter, or normalized frequency for digital filter.
|
||||||
|
@ -1125,11 +1125,11 @@ Else \ref ERROR_CODE_GROUP "code error".
|
||||||
\fn int DSPL_API phase_delay(double* b, double* a, int ord, int flag,
|
\fn int DSPL_API phase_delay(double* b, double* a, int ord, int flag,
|
||||||
double* w, int n, double* tau)
|
double* w, int n, double* tau)
|
||||||
\brief
|
\brief
|
||||||
Расчет задержки цифрового или аналогового фильтра.
|
Расчет фазовой задержки цифрового или аналогового фильтра.
|
||||||
|
|
||||||
Фазовая задержка определяется как:
|
Фазовая задержка определяется как:
|
||||||
\f[
|
\f[
|
||||||
\tau_g(\omega) = - \frac{\Phi(\omega)}{\omega},
|
\tau_{\varphi}(\omega) = - \frac{\Phi(\omega)}{\omega},
|
||||||
\f]
|
\f]
|
||||||
где \f$\Phi(\omega)\f$ -- ФЧХ фильтра, \f$\omega\f$ циктическая частот в случае
|
где \f$\Phi(\omega)\f$ -- ФЧХ фильтра, \f$\omega\f$ циктическая частот в случае
|
||||||
аналогового фильтра, или нормированная частота цифрового фильтра.
|
аналогового фильтра, или нормированная частота цифрового фильтра.
|
||||||
|
|
|
@ -870,8 +870,14 @@ int DSPL_API writetxt_3d(double* x, int nx, double* y, int ny,
|
||||||
for(k = 0; k < ny; k++)
|
for(k = 0; k < ny; k++)
|
||||||
{
|
{
|
||||||
for(n = 0; n < nx; n++)
|
for(n = 0; n < nx; n++)
|
||||||
fprintf(pFile, "%+.12E\t%+.12E\t%+.12E\n",
|
{
|
||||||
x[n], y[k], z[n+k*nx]);
|
if(!isnan(z[n+k*nx]))
|
||||||
|
{
|
||||||
|
fprintf(pFile, "%+.12E\t%+.12E\t%+.12E\n",
|
||||||
|
x[n], y[k], z[n+k*nx]);
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
fprintf(pFile, "\n");
|
fprintf(pFile, "\n");
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -227,6 +227,68 @@ exit_label:
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
|
int DSPL_API mean(double* x, int n, double* m)
|
||||||
|
{
|
||||||
|
int k;
|
||||||
|
double sum;
|
||||||
|
if(!x || !m)
|
||||||
|
return ERROR_PTR;
|
||||||
|
if(n<1)
|
||||||
|
return ERROR_SIZE;
|
||||||
|
|
||||||
|
sum = x[0];
|
||||||
|
for(k = 1; k < n; k++)
|
||||||
|
sum += x[k];
|
||||||
|
|
||||||
|
(*m) = sum / (double)n;
|
||||||
|
return RES_OK;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
|
int DSPL_API mean_cmplx(complex_t* x, int n, complex_t* m)
|
||||||
|
{
|
||||||
|
int k;
|
||||||
|
complex_t sum;
|
||||||
|
if(!x || !m)
|
||||||
|
return ERROR_PTR;
|
||||||
|
if(n<1)
|
||||||
|
return ERROR_SIZE;
|
||||||
|
|
||||||
|
RE(sum) = RE(x[0]);
|
||||||
|
IM(sum) = IM(x[0]);
|
||||||
|
for(k = 1; k < n; k++)
|
||||||
|
{
|
||||||
|
RE(sum) += RE(x[k]);
|
||||||
|
IM(sum) += IM(x[k]);
|
||||||
|
}
|
||||||
|
RE(sum) /= (double)n;
|
||||||
|
IM(sum) /= (double)n;
|
||||||
|
memcpy(m, sum, sizeof(complex_t));
|
||||||
|
|
||||||
|
return RES_OK;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#ifdef DOXYGEN_ENGLISH
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
@ -260,5 +322,60 @@ int DSPL_API minmax(double* x, int n, double* xmin, double* xmax)
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
|
int DSPL_API std(double* x, int n, double* s)
|
||||||
|
{
|
||||||
|
int k, err;
|
||||||
|
double sum, m;
|
||||||
|
err = mean(x, n, &m);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
sum = (x[0] - m) * (x[0] - m);
|
||||||
|
for(k = 1; k < n; k++)
|
||||||
|
sum += (x[k] - m) * (x[k] - m);
|
||||||
|
(*s) = sqrt(sum / (double)(n-1));
|
||||||
|
exit_label:
|
||||||
|
return err;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
|
int DSPL_API std_cmplx(complex_t* x, int n, double* s)
|
||||||
|
{
|
||||||
|
int k, err;
|
||||||
|
complex_t tmp, m;
|
||||||
|
double sum;
|
||||||
|
err = mean_cmplx(x, n, &m);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
|
||||||
|
RE(tmp) = RE(x[0]) - RE(m);
|
||||||
|
IM(tmp) = IM(x[0]) - IM(m);
|
||||||
|
sum = ABSSQR(tmp);
|
||||||
|
for(k = 1; k < n; k++)
|
||||||
|
{
|
||||||
|
RE(tmp) = RE(x[k]) - RE(m);
|
||||||
|
IM(tmp) = IM(x[k]) - IM(m);
|
||||||
|
sum += ABSSQR(tmp);
|
||||||
|
}
|
||||||
|
*s = sqrt(sum / (double)(n-1));
|
||||||
|
exit_label:
|
||||||
|
return err;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -0,0 +1,496 @@
|
||||||
|
/*
|
||||||
|
* 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"
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup ARRAY_GROUP
|
||||||
|
\fn int verif(double* x, double* y, size_t n, double eps, double* err)
|
||||||
|
\brief Real arrays verification
|
||||||
|
|
||||||
|
Function calculates a maximum relative error between two real arrays `x`
|
||||||
|
and `y` (both length equals `n`):
|
||||||
|
|
||||||
|
\f[
|
||||||
|
e = \max \left( \frac{|x(k) - y(k)| }{ |x(k)|} \right),
|
||||||
|
\quad if \quad |x(k)| > 0,
|
||||||
|
\f]
|
||||||
|
or
|
||||||
|
\f[
|
||||||
|
e = \max(|x(k) - y(k)| ), ~\qquad if \quad~|x(k)| = 0,
|
||||||
|
\f]
|
||||||
|
|
||||||
|
This function can be used for algorithms verification if vector `x` is user
|
||||||
|
algorithm result and vector `y` -- reference vector.
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Pointer to the first vector `x`. \n
|
||||||
|
Vector size is `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Pointer to the second vector `y`. \n
|
||||||
|
Vector size is `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] n
|
||||||
|
Size of vectors `x` and `y`. \n \n
|
||||||
|
|
||||||
|
\param[in] eps
|
||||||
|
Relative error threshold. \n
|
||||||
|
If error less than `eps`, then function returns
|
||||||
|
`DSPL_VERIF_SUCCESS`, else `DSPL_VERIF_FAILED`. \n \n
|
||||||
|
|
||||||
|
\param[in, out] err
|
||||||
|
Pointer to the variable which keep
|
||||||
|
maximum relative error. \n
|
||||||
|
Pointer can be `NULL`, maximum error will not be returned
|
||||||
|
in this case. \n \n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`DSPL_VERIF_SUCCESS` if maximum relative error less than `eps`. \n
|
||||||
|
Otherwise `DSPL_VERIF_FAILED`.
|
||||||
|
|
||||||
|
\author Sergey Bakhurin www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup ARRAY_GROUP
|
||||||
|
\fn int verif(double* x, double* y, size_t n, double eps, double* err)
|
||||||
|
\brief Верификация вещественных массивов
|
||||||
|
|
||||||
|
Функция производит расчет максимальной относительной ошибки между вещественными
|
||||||
|
векторами `x` и `y` равной длины `n`:
|
||||||
|
|
||||||
|
\f[
|
||||||
|
e = \max \left( \frac{|x(k) - y(k)| }{ |x(k)|} \right), \quad \quad |x(k)| > 0,
|
||||||
|
\f]
|
||||||
|
или
|
||||||
|
\f[
|
||||||
|
e = \max(|x(k) - y(k)| ), \qquad \quad~|x(k)| = 0,
|
||||||
|
\f]
|
||||||
|
и возвращает `DSPL_VERIF_SUCCESS` если
|
||||||
|
разница \f$ e\f$ меньше `eps`.
|
||||||
|
В противном случае возвращает `DSPL_VERIF_FAILED`. \n
|
||||||
|
Данная функция используется для верификации работы алгоритмов если вектор `x`
|
||||||
|
результат работы алгоритма пользователя, а `y` -- результат работы этого же
|
||||||
|
алгоритма сторонней функцией.
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Указатель на первый вектор `x`. \n
|
||||||
|
Размер вектора `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Указатель на второй вектор `y`. \n
|
||||||
|
Размер вектора `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] n
|
||||||
|
Размер векторов `x` и `y`. \n \n
|
||||||
|
|
||||||
|
\param[in] eps
|
||||||
|
Допустимая относительная ошибка. \n
|
||||||
|
Если максимальная относительная ошибка меньше `eps`, то функция возвращает
|
||||||
|
`DSPL_VERIF_SUCCESS`, в противном случае `DSPL_VERIF_FAILED`. \n \n
|
||||||
|
|
||||||
|
\param[in, out] err
|
||||||
|
Указатель на переменную максимальной относительной ошибки. \n
|
||||||
|
По данному адресу будет записано значение максимальной относительной ошибки. \n
|
||||||
|
Указатель может быть `NULL`, значение ошибки в этом случае
|
||||||
|
не возвращается. \n \n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`DSPL_VERIF_SUCCESS` если относительная ошибка меньше `eps`. \n
|
||||||
|
В противном случае `DSPL_VERIF_FAILED`.
|
||||||
|
|
||||||
|
\author Бахурин Сергей www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
|
#endif
|
||||||
|
int DSPL_API verif(double* x, double* y, size_t n, double eps, double* err)
|
||||||
|
{
|
||||||
|
double d, maxd;
|
||||||
|
size_t k;
|
||||||
|
int res;
|
||||||
|
if(!x || !y)
|
||||||
|
return ERROR_PTR;
|
||||||
|
if(n < 1)
|
||||||
|
return ERROR_SIZE;
|
||||||
|
if(eps <= 0.0 )
|
||||||
|
return ERROR_NEGATIVE;
|
||||||
|
|
||||||
|
maxd = -100.0;
|
||||||
|
|
||||||
|
for(k = 0; k < n; k++)
|
||||||
|
{
|
||||||
|
d = fabs(x[k] - y[k]);
|
||||||
|
if(fabs(x[k]) > 0.0)
|
||||||
|
{
|
||||||
|
d = d / fabs(x[k]);
|
||||||
|
if(d > maxd)
|
||||||
|
maxd = d;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if(err)
|
||||||
|
*err = maxd;
|
||||||
|
|
||||||
|
if(maxd > eps)
|
||||||
|
res = DSPL_VERIF_FAILED;
|
||||||
|
else
|
||||||
|
res = DSPL_VERIF_SUCCESS;
|
||||||
|
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
|
int DSPL_API verif_data_gen(int len, int type, char* fn)
|
||||||
|
{
|
||||||
|
double *pd = NULL;
|
||||||
|
complex_t *pc = NULL;
|
||||||
|
random_t rnd = {0};
|
||||||
|
int err;
|
||||||
|
if(len < 1)
|
||||||
|
return ERROR_SIZE;
|
||||||
|
if(!fn)
|
||||||
|
return ERROR_FNAME;
|
||||||
|
|
||||||
|
err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
switch(type & DAT_MASK)
|
||||||
|
{
|
||||||
|
case DAT_DOUBLE:
|
||||||
|
pd = (double*)malloc(len*sizeof(double));
|
||||||
|
if(!pd)
|
||||||
|
{
|
||||||
|
err = ERROR_MALLOC;
|
||||||
|
goto exit_label;
|
||||||
|
}
|
||||||
|
err = randn(pd, len, 1.0, 10.0, &rnd);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
|
||||||
|
err = writebin(pd, len, type, fn);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
break;
|
||||||
|
case DAT_COMPLEX:
|
||||||
|
pc = (complex_t*)malloc(len*sizeof(complex_t));
|
||||||
|
if(!pc)
|
||||||
|
{
|
||||||
|
err = ERROR_MALLOC;
|
||||||
|
goto exit_label;
|
||||||
|
}
|
||||||
|
err = randn((double*)pc, 2*len, 1.0, 10.0, &rnd);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
|
||||||
|
err = writebin(pc, len, type, fn);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
err = ERROR_DAT_TYPE;
|
||||||
|
}
|
||||||
|
|
||||||
|
exit_label:
|
||||||
|
if(pd)
|
||||||
|
free(pd);
|
||||||
|
if(pc)
|
||||||
|
free(pc);
|
||||||
|
return err;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup ARRAY_GROUP
|
||||||
|
\fn int verif_cmplx(complex_t* x, complex_t* y, size_t n,
|
||||||
|
double eps, double* err)
|
||||||
|
\brief
|
||||||
|
Complex arrays verification
|
||||||
|
|
||||||
|
Function calculates a maximum relative error between two complex arrays `x`
|
||||||
|
and `y` (both length equals `n`):
|
||||||
|
|
||||||
|
\f[
|
||||||
|
e = \max \left( \frac{|x(k) - y(k)| }{ |x(k)|} \right),
|
||||||
|
\quad if \quad |x(k)| > 0,
|
||||||
|
\f]
|
||||||
|
or
|
||||||
|
\f[
|
||||||
|
e = \max(|x(k) - y(k)| ), ~\qquad if \quad~|x(k)| = 0,
|
||||||
|
\f]
|
||||||
|
Return `DSPL_VERIF_SUCCESS` if maximum relative error \f$ e\f$ less than `eps`.
|
||||||
|
Else returns `DSPL_VERIF_FAILED`. \n
|
||||||
|
|
||||||
|
This function can be used for algorithms verification if vector `x` is user
|
||||||
|
algorithm result and vector `y` -- reference vector.
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Pointer to the first vector `x`. \n
|
||||||
|
Vector size is `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Pointer to the second vector `y`. \n
|
||||||
|
Vector size is `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] n
|
||||||
|
Size of vectors `x` and `y`. \n \n
|
||||||
|
|
||||||
|
\param[in] eps
|
||||||
|
Relative error threshold. \n
|
||||||
|
If error less than `eps`, then function returns
|
||||||
|
`DSPL_VERIF_SUCCESS`, else `DSPL_VERIF_FAILED`. \n \n
|
||||||
|
|
||||||
|
\param[in, out] err
|
||||||
|
Pointer to the variable which keep
|
||||||
|
maximum relative error. \n
|
||||||
|
Pointer can be `NULL`, maximum error will not be returned
|
||||||
|
in this case. \n \n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`DSPL_VERIF_SUCCESS` if maximum relative error less than `eps`. \n
|
||||||
|
Otherwise `DSPL_VERIF_FAILED`.
|
||||||
|
|
||||||
|
\author
|
||||||
|
Sergey Bakhurin
|
||||||
|
www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup ARRAY_GROUP
|
||||||
|
\fn int verif_cmplx(complex_t* x, complex_t* y, size_t n,
|
||||||
|
double eps, double* err)
|
||||||
|
\brief Верификация комплексных массивов
|
||||||
|
|
||||||
|
Функция производит расчет максимальной относительной ошибки между комплексными
|
||||||
|
векторами `x` и `y` равной длины `n`:
|
||||||
|
|
||||||
|
\f[
|
||||||
|
e = \max \left( \frac{|x(k) - y(k)|}{|x(k)|} \right), \quad \quad |x(k)| > 0,
|
||||||
|
\f]
|
||||||
|
или
|
||||||
|
\f[
|
||||||
|
e = \max(|x(k) - y(k)| ), ~\qquad \quad~|x(k)| = 0,
|
||||||
|
\f]
|
||||||
|
и возвращает `DSPL_VERIF_SUCCESS` если
|
||||||
|
разница \f$ e\f$ меньше `eps`.
|
||||||
|
В противном случае возвращает `DSPL_VERIF_FAILED`. \n
|
||||||
|
Данная функция используется для верификации работы алгоритмов если вектор `x`
|
||||||
|
результат работы алгоритма пользователя, а `y` -- результат работы этого же
|
||||||
|
алгоритма сторонней функцией.
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Указатель на первый вектор `x`. \n
|
||||||
|
Размер вектора `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Указатель на второй вектор `y`. \n
|
||||||
|
Размер вектора `[n x 1]`. \n \n
|
||||||
|
|
||||||
|
\param[in] n
|
||||||
|
Размер векторов `x` и `y`. \n \n
|
||||||
|
|
||||||
|
\param[in] eps
|
||||||
|
Допустимая относительная ошибка. \n
|
||||||
|
Если максимальная относительная ошибка меньше `eps`, то функция возвращает
|
||||||
|
`DSPL_VERIF_SUCCESS`, в противном случае `DSPL_VERIF_FAILED`. \n \n
|
||||||
|
|
||||||
|
\param[in, out] err
|
||||||
|
Указатель на переменную максимальной относительной ошибки. \n
|
||||||
|
По данному адресу будет записано значение максимальной относительной ошибки. \n
|
||||||
|
Указатель может быть `NULL`, значение ошибки в этом
|
||||||
|
случае не возвращается. \n \n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`DSPL_VERIF_SUCCESS` если функция выполнена успешно. \n
|
||||||
|
В противном случае `DSPL_VERIF_FAILED`.
|
||||||
|
|
||||||
|
\author
|
||||||
|
Бахурин Сергей www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
|
#endif
|
||||||
|
int DSPL_API verif_cmplx(complex_t* x, complex_t* y, size_t n,
|
||||||
|
double eps, double* err)
|
||||||
|
{
|
||||||
|
|
||||||
|
complex_t d;
|
||||||
|
double mx, md, maxd;
|
||||||
|
size_t k;
|
||||||
|
int res;
|
||||||
|
if(!x || !y)
|
||||||
|
return ERROR_PTR;
|
||||||
|
if(n < 1)
|
||||||
|
return ERROR_SIZE;
|
||||||
|
if(eps <= 0.0 )
|
||||||
|
return ERROR_NEGATIVE;
|
||||||
|
|
||||||
|
maxd = -100.0;
|
||||||
|
|
||||||
|
for(k = 0; k < n; k++)
|
||||||
|
{
|
||||||
|
RE(d) = RE(x[k]) - RE(y[k]);
|
||||||
|
IM(d) = IM(x[k]) - IM(y[k]);
|
||||||
|
md = ABS(d);
|
||||||
|
mx = ABS(x[k]);
|
||||||
|
if(mx > 0.0)
|
||||||
|
{
|
||||||
|
md = md / mx;
|
||||||
|
if(md > maxd)
|
||||||
|
maxd = md;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if(err)
|
||||||
|
*err = maxd;
|
||||||
|
|
||||||
|
if(maxd > eps)
|
||||||
|
res = DSPL_VERIF_FAILED;
|
||||||
|
else
|
||||||
|
res = DSPL_VERIF_SUCCESS;
|
||||||
|
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
|
void DSPL_API verif_str(double* yout, int nout,
|
||||||
|
char* str_msg, char* outfn, char* logfn)
|
||||||
|
{
|
||||||
|
char str[VERIF_STR_BUF] = {0};
|
||||||
|
char msg[VERIF_STR_BUF] = {0};
|
||||||
|
double *y = NULL;
|
||||||
|
double derr = 0.0;
|
||||||
|
int n, verr, type;
|
||||||
|
|
||||||
|
sprintf(str, "%s", str_msg);
|
||||||
|
while(strlen(str) < VERIF_STR_LEN)
|
||||||
|
str[strlen(str)] = VERIF_CHAR_POINT;
|
||||||
|
|
||||||
|
readbin(outfn, (void**)(&y), &n, &type);
|
||||||
|
|
||||||
|
if(nout != n)
|
||||||
|
{
|
||||||
|
sprintf(msg, "FAILED (out size [%d] neq [%d])", n, nout);
|
||||||
|
strcat(str, msg);
|
||||||
|
addlog(str, logfn);
|
||||||
|
printf("%s\n", str);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(type!=DAT_DOUBLE)
|
||||||
|
{
|
||||||
|
sprintf(msg, "FAILED (type is complex)");
|
||||||
|
strcat(str, msg);
|
||||||
|
addlog(str, logfn);
|
||||||
|
printf("%s\n", str);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
verr = verif(yout, y, nout, VERIF_LEVEL_DOUBLE, &derr);
|
||||||
|
if(verr == DSPL_VERIF_SUCCESS)
|
||||||
|
sprintf(msg, "ok (err = %12.4E)", derr);
|
||||||
|
else
|
||||||
|
sprintf(msg, "FAILED (err = %12.4E)", derr);
|
||||||
|
strcat(str, msg);
|
||||||
|
addlog(str, logfn);
|
||||||
|
printf("%s\n", str);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
|
void DSPL_API verif_str_cmplx(complex_t* yout, int nout,
|
||||||
|
char* str_msg, char* outfn, char* logfn)
|
||||||
|
{
|
||||||
|
char str[VERIF_STR_BUF] = {0};
|
||||||
|
char msg[VERIF_STR_BUF] = {0};
|
||||||
|
complex_t *y = NULL;
|
||||||
|
double derr = 0.0;
|
||||||
|
int n, verr, type;
|
||||||
|
|
||||||
|
sprintf(str, "%s", str_msg);
|
||||||
|
while(strlen(str) < VERIF_STR_LEN)
|
||||||
|
str[strlen(str)] = VERIF_CHAR_POINT;
|
||||||
|
|
||||||
|
readbin(outfn, (void**)(&y), &n, &type);
|
||||||
|
|
||||||
|
if(nout != n)
|
||||||
|
{
|
||||||
|
sprintf(msg, "FAILED (out size [%d] neq [%d])", n, nout);
|
||||||
|
strcat(str, msg);
|
||||||
|
addlog(str, logfn);
|
||||||
|
printf("%s\n", str);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(type!=DAT_COMPLEX)
|
||||||
|
{
|
||||||
|
sprintf(msg, "FAILED (type is complex)");
|
||||||
|
strcat(str, msg);
|
||||||
|
addlog(str, logfn);
|
||||||
|
printf("%s\n", str);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
verr = verif_cmplx(yout, y, nout, VERIF_LEVEL_DOUBLE, &derr);
|
||||||
|
if(verr == DSPL_VERIF_SUCCESS)
|
||||||
|
sprintf(msg, "ok (err = %12.4E)", derr);
|
||||||
|
else
|
||||||
|
sprintf(msg, "FAILED (err = %12.4E)", derr);
|
||||||
|
strcat(str, msg);
|
||||||
|
addlog(str, logfn);
|
||||||
|
printf("%s\n", str);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
362
dspl/src/xcorr.c
362
dspl/src/xcorr.c
|
@ -36,29 +36,203 @@ int xcorr_scale_cmplx(complex_t* x, int nd, int flag);
|
||||||
|
|
||||||
|
|
||||||
#ifdef DOXYGEN_ENGLISH
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup SPEC_MATH_STAT_GROUP
|
||||||
|
\fn int xcorr(double* x, int nx, double* y, int ny,
|
||||||
|
int flag, int nr, double* r, double* t)
|
||||||
|
\brief Estimates the cross-correlation vector for real
|
||||||
|
discrete-time sequences `x` and `y`.
|
||||||
|
|
||||||
|
Estimate the cross correlation \f$\widehat{r}_{xy}(k)\f$ of vector arguments
|
||||||
|
`x` and `y` or estimate autocorrelation vector if \f$x = y \f$.
|
||||||
|
|
||||||
|
Cross-correlation vector is defined as:
|
||||||
|
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N-k} \sum\limits_{n = 0}^{N-k-1} x(n+k)y^*(n), \qquad 0 \leq k <N;
|
||||||
|
\f]
|
||||||
|
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N+k} \sum\limits_{n = 0}^{N+k-1} y^*(n-k)x(n), \qquad -N < k < 0.
|
||||||
|
\f]
|
||||||
|
here \f$ N = \max(n_x, n_y) \f$.
|
||||||
|
|
||||||
|
This function uses the FFT algorithm to estimate the scaled correlation vector:
|
||||||
|
\f[
|
||||||
|
\breve{r}_{xy}(k) = \operatorname{IFFT}\Big[
|
||||||
|
\operatorname{FFT}\big[ \mathbf{x} \big]
|
||||||
|
\operatorname{FFT}^*\big[ \mathbf{y} \big]
|
||||||
|
\Big]
|
||||||
|
\f]
|
||||||
|
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Pointer to the discrete-time vector `x`. \n
|
||||||
|
Vector size is `[nx x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nx
|
||||||
|
Size of vector `x`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Pointer to the discrete-time vector `y`. \n
|
||||||
|
Vector size is `[ny x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] ny
|
||||||
|
Size of vector `y`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] flag
|
||||||
|
Flag specifies the type of scaling applied to the correlation vector
|
||||||
|
\f$\breve{r}_{xy}(k)\f$.\n
|
||||||
|
Is one of:\n
|
||||||
|
`DSPL_XCORR_NOSCALE` unscaled correlation vector from IFFT output
|
||||||
|
\f$\breve{r}_{xy}(k)\f$;\n
|
||||||
|
`DSPL_XCORR_BIASED` biased correlation vector \f$\breve{r}_{xy}(k)/N \f$;\n
|
||||||
|
`DSPL_XCORR_UNBIASED` unbiased correlation vector
|
||||||
|
\f$\widehat{r}_{xy}(k) = \frac{\breve{r}_{xy}(k)}{N-|k|} \f$;\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nr
|
||||||
|
Maximum correlation lag.\n
|
||||||
|
Correlation vector \f$\widehat{r}_{xy}(k)\f$ is calculated for
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldots n_r\f$.\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] r
|
||||||
|
Pointer to the cross-correlation or autocorrelation vector. \n
|
||||||
|
Vector size is `[(2*nr+1) x 1]`. \n
|
||||||
|
Memory must be allocated. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] r
|
||||||
|
Pointer to the cross-correlation argument vector
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldors n_r\f$.\n
|
||||||
|
Vector size is `[(2*nr+1) x 1]`. \n
|
||||||
|
Pointer can be `NULL`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`RES_OK` if function returns successfully. \n
|
||||||
|
Else \ref ERROR_CODE_GROUP "code error".
|
||||||
|
|
||||||
|
\author Sergey Bakhurin www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
#endif
|
#endif
|
||||||
#ifdef DOXYGEN_RUSSIAN
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup SPEC_MATH_STAT_GROUP
|
||||||
|
\fn int xcorr(double* x, int nx, double* y, int ny,
|
||||||
|
int flag, int nr, double* r, double* t)
|
||||||
|
\brief Оценка вектора взаимной корреляции для дискретных
|
||||||
|
вещественных последовательностей `x` и `y`.
|
||||||
|
|
||||||
|
Функция производит оценку вектора взаимной корреляции \f$\widehat{r}_{xy}(k)\f$
|
||||||
|
для векторов `x` и `y` или вектора автокорреляции
|
||||||
|
\f$\widehat{r}_{xx}(k)\f$ если \f$x = y \f$.
|
||||||
|
|
||||||
|
Несмещенная оценка вектора взаимной корреляции:
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N-k} \sum\limits_{n = 0}^{N-k-1} x(n+k)y^*(n), \qquad 0 \leq k <N;
|
||||||
|
\f]
|
||||||
|
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N+k} \sum\limits_{n = 0}^{N+k-1} y^*(n-k)x(n), \qquad -N < k < 0.
|
||||||
|
\f]
|
||||||
|
где \f$ N = \max(n_x, n_y) \f$.
|
||||||
|
|
||||||
|
Данная функция использует алгоритм FFT для вычислительной эффективности оценки:
|
||||||
|
\f[
|
||||||
|
\breve{r}_{xy}(k) = \operatorname{IFFT}\Big[
|
||||||
|
\operatorname{FFT}\big[ \mathbf{x} \big]
|
||||||
|
\operatorname{FFT}^*\big[ \mathbf{y} \big]
|
||||||
|
\Big]
|
||||||
|
\f]
|
||||||
|
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Указатель на первую дискретную последовательность `x`. \n
|
||||||
|
Размер вектора `[nx x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nx
|
||||||
|
Размер вектора первой дискретной последовательности `x`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Указатель на вторую дискретную последовательность `y`. \n
|
||||||
|
Размер вектора `[ny x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] ny
|
||||||
|
Размер вектора второй дискретной последовательности `y`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] flag
|
||||||
|
Флаг задает способ масштабирования выходного корреляционного вектора
|
||||||
|
\f$\breve{r}_{xy}(k)\f$.\n
|
||||||
|
Может принимать одно из следующих значений:\n
|
||||||
|
`DSPL_XCORR_NOSCALE` немасштабированный выход алгоритма IFFT
|
||||||
|
\f$\breve{r}_{xy}(k)\f$;\n
|
||||||
|
`DSPL_XCORR_BIASED` Смещенная оценка \f$\breve{r}_{xy}(k)/N \f$;\n
|
||||||
|
`DSPL_XCORR_UNBIASED` Несмещенная оценка
|
||||||
|
\f$\widehat{r}_{xy}(k) = \frac{\breve{r}_{xy}(k)}{N-|k|} \f$;\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nr
|
||||||
|
Диапазон оценки вектора корреляции относительно нуля.\n
|
||||||
|
Вектор \f$\widehat{r}_{xy}(k)\f$ рассчитывается для значений аргумента
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldots n_r\f$.\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] r
|
||||||
|
Указатель на масштабированный вектор взаимной корреляции. \n
|
||||||
|
Размер вектора `[(2*nr+1) x 1]`. \n
|
||||||
|
Память должна быть выделена. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] t
|
||||||
|
Указатель на значения аргумента вектора взаимной корреляции
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldors n_r\f$.\n
|
||||||
|
Размер вектора `[(2*nr+1) x 1]`. \n
|
||||||
|
Указатель может быть `NULL`. В этом случае значения аргумента не возвращаются.\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`RES_OK` Если функция рассчитана успешно. \n
|
||||||
|
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
|
||||||
|
|
||||||
|
\author Бахурин Сергей www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
#endif
|
#endif
|
||||||
int DSPL_API xcorr(double* x, int nx, double* y, int ny,
|
int DSPL_API xcorr(double* x, int nx, double* y, int ny,
|
||||||
int flag, int nr, double* r, double* t)
|
int flag, int nr, double* r, double* t)
|
||||||
{
|
{
|
||||||
fft_t fft = {0};
|
fft_t fft = {0};
|
||||||
int err;
|
int err;
|
||||||
complex_t *cr = (complex_t*)malloc((2 * nr + 1) * sizeof(complex_t));
|
complex_t *cx = NULL;
|
||||||
|
complex_t *cy = NULL;
|
||||||
|
complex_t *cr = NULL;
|
||||||
|
|
||||||
|
cr = (complex_t*)malloc((2 * nr + 1) * sizeof(complex_t));
|
||||||
if(!cr)
|
if(!cr)
|
||||||
{
|
{
|
||||||
err = ERROR_MALLOC;
|
err = ERROR_MALLOC;
|
||||||
goto exit_label;
|
goto exit_label;
|
||||||
}
|
}
|
||||||
complex_t *cx = (complex_t*)malloc( nx * sizeof(complex_t));
|
cx = (complex_t*)malloc( nx * sizeof(complex_t));
|
||||||
if(!cx)
|
if(!cx)
|
||||||
{
|
{
|
||||||
err = ERROR_MALLOC;
|
err = ERROR_MALLOC;
|
||||||
goto exit_label;
|
goto exit_label;
|
||||||
}
|
}
|
||||||
complex_t *cy = (complex_t*)malloc( ny * sizeof(complex_t));
|
cy = (complex_t*)malloc( ny * sizeof(complex_t));
|
||||||
if(!cy)
|
if(!cy)
|
||||||
{
|
{
|
||||||
err = ERROR_MALLOC;
|
err = ERROR_MALLOC;
|
||||||
|
@ -97,10 +271,180 @@ exit_label:
|
||||||
|
|
||||||
|
|
||||||
#ifdef DOXYGEN_ENGLISH
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup SPEC_MATH_STAT_GROUP
|
||||||
|
int xcorr_cmplx(complex_t* x, int nx, complex_t* y, int ny,
|
||||||
|
int flag, int nr, complex_t* r, double* t)
|
||||||
|
\brief Estimates the cross-correlation vector for complex
|
||||||
|
discrete-time sequences `x` and `y`.
|
||||||
|
|
||||||
|
Estimate the cross correlation \f$\widehat{r}_{xy}(k)\f$ of vector arguments
|
||||||
|
`x` and `y` or estimate autocorrelation vector if \f$x = y \f$.
|
||||||
|
|
||||||
|
Cross-correlation vector is defined as:
|
||||||
|
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N-k} \sum\limits_{n = 0}^{N-k-1} x(n+k)y^*(n), \qquad 0 \leq k <N;
|
||||||
|
\f]
|
||||||
|
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N+k} \sum\limits_{n = 0}^{N+k-1} y^*(n-k)x(n), \qquad -N < k < 0.
|
||||||
|
\f]
|
||||||
|
here \f$ N = \max(n_x, n_y) \f$.
|
||||||
|
|
||||||
|
This function uses the FFT algorithm to estimate the scaled correlation vector:
|
||||||
|
\f[
|
||||||
|
\breve{r}_{xy}(k) = \operatorname{IFFT}\Big[
|
||||||
|
\operatorname{FFT}\big[ \mathbf{x} \big]
|
||||||
|
\operatorname{FFT}^*\big[ \mathbf{y} \big]
|
||||||
|
\Big]
|
||||||
|
\f]
|
||||||
|
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Pointer to the discrete-time vector `x`. \n
|
||||||
|
Vector size is `[nx x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nx
|
||||||
|
Size of vector `x`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Pointer to the discrete-time vector `y`. \n
|
||||||
|
Vector size is `[ny x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] ny
|
||||||
|
Size of vector `y`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] flag
|
||||||
|
Flag specifies the type of scaling applied to the correlation vector
|
||||||
|
\f$\breve{r}_{xy}(k)\f$.\n
|
||||||
|
Is one of:\n
|
||||||
|
`DSPL_XCORR_NOSCALE` unscaled correlation vector from IFFT output
|
||||||
|
\f$\breve{r}_{xy}(k)\f$;\n
|
||||||
|
`DSPL_XCORR_BIASED` biased correlation vector \f$\breve{r}_{xy}(k)/N \f$;\n
|
||||||
|
`DSPL_XCORR_UNBIASED` unbiased correlation vector
|
||||||
|
\f$\widehat{r}_{xy}(k) = \frac{\breve{r}_{xy}(k)}{N-|k|} \f$;\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nr
|
||||||
|
Maximum correlation lag.\n
|
||||||
|
Correlation vector \f$\widehat{r}_{xy}(k)\f$ is calculated for
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldots n_r\f$.\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] r
|
||||||
|
Pointer to the cross-correlation or autocorrelation vector. \n
|
||||||
|
Vector size is `[(2*nr+1) x 1]`. \n
|
||||||
|
Memory must be allocated. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] r
|
||||||
|
Pointer to the cross-correlation argument vector
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldors n_r\f$.\n
|
||||||
|
Vector size is `[(2*nr+1) x 1]`. \n
|
||||||
|
Pointer can be `NULL`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`RES_OK` if function returns successfully. \n
|
||||||
|
Else \ref ERROR_CODE_GROUP "code error".
|
||||||
|
|
||||||
|
\author Sergey Bakhurin www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
#endif
|
#endif
|
||||||
#ifdef DOXYGEN_RUSSIAN
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
/*! ****************************************************************************
|
||||||
|
\ingroup SPEC_MATH_STAT_GROUP
|
||||||
|
int xcorr_cmplx(complex_t* x, int nx, complex_t* y, int ny,
|
||||||
|
int flag, int nr, complex_t* r, double* t)
|
||||||
|
\brief Оценка вектора взаимной корреляции для дискретных
|
||||||
|
комплексных последовательностей `x` и `y`.
|
||||||
|
|
||||||
|
Функция производит оценку вектора взаимной корреляции \f$\widehat{r}_{xy}(k)\f$
|
||||||
|
для векторов `x` и `y` или вектора автокорреляции
|
||||||
|
\f$\widehat{r}_{xx}(k)\f$ если \f$x = y \f$.
|
||||||
|
|
||||||
|
Несмещенная оценка вектора взаимной корреляции:
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N-k} \sum\limits_{n = 0}^{N-k-1} x(n+k)y^*(n), \qquad 0 \leq k <N;
|
||||||
|
\f]
|
||||||
|
|
||||||
|
\f[
|
||||||
|
\widehat{r}_{xy}(k) =
|
||||||
|
\frac{1}{N+k} \sum\limits_{n = 0}^{N+k-1} y^*(n-k)x(n), \qquad -N < k < 0.
|
||||||
|
\f]
|
||||||
|
где \f$ N = \max(n_x, n_y) \f$.
|
||||||
|
|
||||||
|
Данная функция использует алгоритм FFT для вычислительной эффективности оценки:
|
||||||
|
\f[
|
||||||
|
\breve{r}_{xy}(k) = \operatorname{IFFT}\Big[
|
||||||
|
\operatorname{FFT}\big[ \mathbf{x} \big]
|
||||||
|
\operatorname{FFT}^*\big[ \mathbf{y} \big]
|
||||||
|
\Big]
|
||||||
|
\f]
|
||||||
|
|
||||||
|
|
||||||
|
\param[in] x
|
||||||
|
Указатель на первую дискретную последовательность `x`. \n
|
||||||
|
Размер вектора `[nx x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nx
|
||||||
|
Размер вектора первой дискретной последовательности `x`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] y
|
||||||
|
Указатель на вторую дискретную последовательность `y`. \n
|
||||||
|
Размер вектора `[ny x 1]`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] ny
|
||||||
|
Размер вектора второй дискретной последовательности `y`. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] flag
|
||||||
|
Флаг задает способ масштабирования выходного корреляционного вектора
|
||||||
|
\f$\breve{r}_{xy}(k)\f$.\n
|
||||||
|
Может принимать одно из следующих значений:\n
|
||||||
|
`DSPL_XCORR_NOSCALE` немасштабированный выход алгоритма IFFT
|
||||||
|
\f$\breve{r}_{xy}(k)\f$;\n
|
||||||
|
`DSPL_XCORR_BIASED` Смещенная оценка \f$\breve{r}_{xy}(k)/N \f$;\n
|
||||||
|
`DSPL_XCORR_UNBIASED` Несмещенная оценка
|
||||||
|
\f$\widehat{r}_{xy}(k) = \frac{\breve{r}_{xy}(k)}{N-|k|} \f$;\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[in] nr
|
||||||
|
Диапазон оценки вектора корреляции относительно нуля.\n
|
||||||
|
Вектор \f$\widehat{r}_{xy}(k)\f$ рассчитывается для значений аргумента
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldots n_r\f$.\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] r
|
||||||
|
Указатель на масштабированный вектор взаимной корреляции. \n
|
||||||
|
Размер вектора `[(2*nr+1) x 1]`. \n
|
||||||
|
Память должна быть выделена. \n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\param[out] t
|
||||||
|
Указатель на значения аргумента вектора взаимной корреляции
|
||||||
|
\f$ k= -n_r,\,\, -n_r +1, \ldors n_r\f$.\n
|
||||||
|
Размер вектора `[(2*nr+1) x 1]`. \n
|
||||||
|
Указатель может быть `NULL`. В этом случае значения аргумента не возвращаются.\n
|
||||||
|
\n
|
||||||
|
|
||||||
|
\return
|
||||||
|
`RES_OK` Если функция рассчитана успешно. \n
|
||||||
|
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
|
||||||
|
|
||||||
|
\author Бахурин Сергей www.dsplib.org
|
||||||
|
***************************************************************************** */
|
||||||
#endif
|
#endif
|
||||||
int DSPL_API xcorr_cmplx(complex_t* x, int nx, complex_t* y, int ny,
|
int DSPL_API xcorr_cmplx(complex_t* x, int nx, complex_t* y, int ny,
|
||||||
int flag, int nr, complex_t* r, double* t)
|
int flag, int nr, complex_t* r, double* t)
|
||||||
|
@ -115,7 +459,12 @@ int DSPL_API xcorr_cmplx(complex_t* x, int nx, complex_t* y, int ny,
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
int xcorr_get_lag_cmplx(complex_t* x, int nd, int nr, complex_t* r, double* t)
|
int xcorr_get_lag_cmplx(complex_t* x, int nd, int nr, complex_t* r, double* t)
|
||||||
{
|
{
|
||||||
int i;
|
int i;
|
||||||
|
@ -148,7 +497,7 @@ int xcorr_get_lag_cmplx(complex_t* x, int nd, int nr, complex_t* r, double* t)
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
int xcorr_krn(complex_t* x, int nx, complex_t* y, int ny, fft_t* pfft,
|
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 flag, int nr, complex_t* r, double* t)
|
||||||
{
|
{
|
||||||
complex_t *px = NULL;
|
complex_t *px = NULL;
|
||||||
complex_t *py = NULL;
|
complex_t *py = NULL;
|
||||||
|
@ -351,7 +700,12 @@ int xcorr_fft_size(int nx, int ny, int* pnfft, int* pndata)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef DOXYGEN_ENGLISH
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#ifdef DOXYGEN_RUSSIAN
|
||||||
|
|
||||||
|
#endif
|
||||||
int xcorr_scale_cmplx(complex_t* x, int nd, int flag)
|
int xcorr_scale_cmplx(complex_t* x, int nd, int flag)
|
||||||
{
|
{
|
||||||
int i;
|
int i;
|
||||||
|
|
|
@ -140,6 +140,8 @@ p_matrix_print_cmplx matrix_print_cmplx ;
|
||||||
p_matrix_transpose matrix_transpose ;
|
p_matrix_transpose matrix_transpose ;
|
||||||
p_matrix_transpose_cmplx matrix_transpose_cmplx ;
|
p_matrix_transpose_cmplx matrix_transpose_cmplx ;
|
||||||
p_matrix_transpose_hermite matrix_transpose_hermite ;
|
p_matrix_transpose_hermite matrix_transpose_hermite ;
|
||||||
|
p_mean mean ;
|
||||||
|
p_mean_cmplx mean_cmplx ;
|
||||||
p_minmax minmax ;
|
p_minmax minmax ;
|
||||||
|
|
||||||
p_ones ones ;
|
p_ones ones ;
|
||||||
|
@ -166,6 +168,8 @@ p_sin_cmplx sin_cmplx ;
|
||||||
p_sinc sinc ;
|
p_sinc sinc ;
|
||||||
p_sine_int sine_int ;
|
p_sine_int sine_int ;
|
||||||
p_sqrt_cmplx sqrt_cmplx ;
|
p_sqrt_cmplx sqrt_cmplx ;
|
||||||
|
p_std std ;
|
||||||
|
p_std_cmplx std_cmplx ;
|
||||||
|
|
||||||
p_trapint trapint ;
|
p_trapint trapint ;
|
||||||
p_trapint_cmplx trapint_cmplx ;
|
p_trapint_cmplx trapint_cmplx ;
|
||||||
|
@ -174,7 +178,10 @@ p_unwrap unwrap ;
|
||||||
|
|
||||||
p_vector_dot vector_dot ;
|
p_vector_dot vector_dot ;
|
||||||
p_verif verif ;
|
p_verif verif ;
|
||||||
|
p_verif_data_gen verif_data_gen ;
|
||||||
p_verif_cmplx verif_cmplx ;
|
p_verif_cmplx verif_cmplx ;
|
||||||
|
p_verif_str verif_str ;
|
||||||
|
p_verif_str_cmplx verif_str_cmplx ;
|
||||||
|
|
||||||
p_window window ;
|
p_window window ;
|
||||||
p_writebin writebin ;
|
p_writebin writebin ;
|
||||||
|
@ -447,6 +454,8 @@ void* dspl_load()
|
||||||
LOAD_FUNC(matrix_transpose);
|
LOAD_FUNC(matrix_transpose);
|
||||||
LOAD_FUNC(matrix_transpose_cmplx);
|
LOAD_FUNC(matrix_transpose_cmplx);
|
||||||
LOAD_FUNC(matrix_transpose_hermite);
|
LOAD_FUNC(matrix_transpose_hermite);
|
||||||
|
LOAD_FUNC(mean);
|
||||||
|
LOAD_FUNC(mean_cmplx);
|
||||||
LOAD_FUNC(minmax);
|
LOAD_FUNC(minmax);
|
||||||
|
|
||||||
LOAD_FUNC(ones);
|
LOAD_FUNC(ones);
|
||||||
|
@ -473,6 +482,8 @@ void* dspl_load()
|
||||||
LOAD_FUNC(sinc);
|
LOAD_FUNC(sinc);
|
||||||
LOAD_FUNC(sine_int);
|
LOAD_FUNC(sine_int);
|
||||||
LOAD_FUNC(sqrt_cmplx);
|
LOAD_FUNC(sqrt_cmplx);
|
||||||
|
LOAD_FUNC(std);
|
||||||
|
LOAD_FUNC(std_cmplx);
|
||||||
|
|
||||||
LOAD_FUNC(trapint);
|
LOAD_FUNC(trapint);
|
||||||
LOAD_FUNC(trapint_cmplx);
|
LOAD_FUNC(trapint_cmplx);
|
||||||
|
@ -481,7 +492,10 @@ void* dspl_load()
|
||||||
|
|
||||||
LOAD_FUNC(vector_dot);
|
LOAD_FUNC(vector_dot);
|
||||||
LOAD_FUNC(verif);
|
LOAD_FUNC(verif);
|
||||||
|
LOAD_FUNC(verif_data_gen);
|
||||||
LOAD_FUNC(verif_cmplx);
|
LOAD_FUNC(verif_cmplx);
|
||||||
|
LOAD_FUNC(verif_str);
|
||||||
|
LOAD_FUNC(verif_str_cmplx);
|
||||||
|
|
||||||
LOAD_FUNC(window);
|
LOAD_FUNC(window);
|
||||||
LOAD_FUNC(writebin);
|
LOAD_FUNC(writebin);
|
||||||
|
|
|
@ -635,6 +635,15 @@ Variable `z = 1-2j`, here `j` - imaginary unit, but variable `y = 5`.
|
||||||
#define PLOT_HOLD 0x00000001
|
#define PLOT_HOLD 0x00000001
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#define VERIF_STR_BUF 128
|
||||||
|
#define VERIF_STR_LEN 48
|
||||||
|
#define VERIF_CHAR_POINT 46
|
||||||
|
#define VERIF_LEVEL_COMPLEX 1E-11
|
||||||
|
#define VERIF_LEVEL_DOUBLE 1E-12
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#ifdef __cplusplus
|
#ifdef __cplusplus
|
||||||
extern "C" {
|
extern "C" {
|
||||||
#endif
|
#endif
|
||||||
|
@ -1214,12 +1223,20 @@ DECLARE_FUNC(int, matrix_transpose_hermite, complex_t* a
|
||||||
COMMA int m
|
COMMA int m
|
||||||
COMMA complex_t* b);
|
COMMA complex_t* b);
|
||||||
/*----------------------------------------------------------------------------*/
|
/*----------------------------------------------------------------------------*/
|
||||||
|
DECLARE_FUNC(int, mean, double* x
|
||||||
|
COMMA int n
|
||||||
|
COMMA double* m);
|
||||||
|
/*----------------------------------------------------------------------------*/
|
||||||
|
DECLARE_FUNC(int, mean_cmplx, complex_t* x
|
||||||
|
COMMA int n
|
||||||
|
COMMA complex_t* m);
|
||||||
|
/*----------------------------------------------------------------------------*/
|
||||||
DECLARE_FUNC(int, minmax, double* x
|
DECLARE_FUNC(int, minmax, double* x
|
||||||
COMMA int n
|
COMMA int n
|
||||||
COMMA double* xmin
|
COMMA double* xmin
|
||||||
COMMA double* xmax);
|
COMMA double* xmax);
|
||||||
/*----------------------------------------------------------------------------*/
|
/*----------------------------------------------------------------------------*/
|
||||||
DECLARE_FUNC(int, ones, double* x
|
DECLARE_FUNC(int, ones, double* x
|
||||||
COMMA int n);
|
COMMA int n);
|
||||||
/*----------------------------------------------------------------------------*/
|
/*----------------------------------------------------------------------------*/
|
||||||
DECLARE_FUNC(int, phase_delay, double* b
|
DECLARE_FUNC(int, phase_delay, double* b
|
||||||
|
@ -1330,6 +1347,14 @@ DECLARE_FUNC(int, sqrt_cmplx, complex_t*
|
||||||
COMMA int
|
COMMA int
|
||||||
COMMA complex_t*);
|
COMMA complex_t*);
|
||||||
/*----------------------------------------------------------------------------*/
|
/*----------------------------------------------------------------------------*/
|
||||||
|
DECLARE_FUNC(int, std, double* x
|
||||||
|
COMMA int n
|
||||||
|
COMMA double* s);
|
||||||
|
/*----------------------------------------------------------------------------*/
|
||||||
|
DECLARE_FUNC(int, std_cmplx, complex_t* x
|
||||||
|
COMMA int n
|
||||||
|
COMMA double* s);
|
||||||
|
/*----------------------------------------------------------------------------*/
|
||||||
DECLARE_FUNC(int, trapint, double*
|
DECLARE_FUNC(int, trapint, double*
|
||||||
COMMA double*
|
COMMA double*
|
||||||
COMMA int
|
COMMA int
|
||||||
|
@ -1356,12 +1381,28 @@ DECLARE_FUNC(int, verif, double* x
|
||||||
COMMA double eps
|
COMMA double eps
|
||||||
COMMA double* err);
|
COMMA double* err);
|
||||||
/*----------------------------------------------------------------------------*/
|
/*----------------------------------------------------------------------------*/
|
||||||
|
DECLARE_FUNC(int, verif_data_gen, int len
|
||||||
|
COMMA int type
|
||||||
|
COMMA char* fn);
|
||||||
|
/*----------------------------------------------------------------------------*/
|
||||||
DECLARE_FUNC(int, verif_cmplx, complex_t* x
|
DECLARE_FUNC(int, verif_cmplx, complex_t* x
|
||||||
COMMA complex_t* y
|
COMMA complex_t* y
|
||||||
COMMA size_t n
|
COMMA size_t n
|
||||||
COMMA double eps
|
COMMA double eps
|
||||||
COMMA double* err);
|
COMMA double* err);
|
||||||
/*----------------------------------------------------------------------------*/
|
/*----------------------------------------------------------------------------*/
|
||||||
|
DECLARE_FUNC(void, verif_str, double* yout
|
||||||
|
COMMA int nout
|
||||||
|
COMMA char* str_msg
|
||||||
|
COMMA char* outfn
|
||||||
|
COMMA char* logfn);
|
||||||
|
/*----------------------------------------------------------------------------*/
|
||||||
|
DECLARE_FUNC(void, verif_str_cmplx, complex_t* yout
|
||||||
|
COMMA int nout
|
||||||
|
COMMA char* str_msg
|
||||||
|
COMMA char* outfn
|
||||||
|
COMMA char* logfn);
|
||||||
|
/*----------------------------------------------------------------------------*/
|
||||||
DECLARE_FUNC(int, window, double* w
|
DECLARE_FUNC(int, window, double* w
|
||||||
COMMA int n
|
COMMA int n
|
||||||
COMMA int win_type
|
COMMA int win_type
|
||||||
|
|
|
@ -0,0 +1,3 @@
|
||||||
|
*.txt
|
||||||
|
*.bin
|
||||||
|
*.dat
|
|
@ -0,0 +1,71 @@
|
||||||
|
clear all; close all; clc;
|
||||||
|
|
||||||
|
|
||||||
|
N = 262144*16;
|
||||||
|
|
||||||
|
x = randn(N,1) + 1i * randn(N,1);
|
||||||
|
|
||||||
|
|
||||||
|
k = 4;
|
||||||
|
s = N;
|
||||||
|
|
||||||
|
for j = 1:21
|
||||||
|
z = x(1:s);
|
||||||
|
t = tic();
|
||||||
|
for i = 1:k
|
||||||
|
y = fft(z);
|
||||||
|
endfor
|
||||||
|
dt = toc(t) / k;
|
||||||
|
mflops(22-j) = 5 * s * log2(s) / dt / 1E6 ;
|
||||||
|
size(22-j) = s;
|
||||||
|
fprintf("size = %8d Mflops = %12.3f\n", s, mflops(22-j));
|
||||||
|
k = k * 1.7;
|
||||||
|
s = s / 2;
|
||||||
|
endfor
|
||||||
|
|
||||||
|
|
||||||
|
dspl = [1204.630392
|
||||||
|
1283.970612
|
||||||
|
1586.347958
|
||||||
|
1707.107097
|
||||||
|
1866.109831
|
||||||
|
1837.307509
|
||||||
|
2366.785829
|
||||||
|
2302.925874
|
||||||
|
2388.456514
|
||||||
|
2113.451546
|
||||||
|
3090.904615
|
||||||
|
2979.596190
|
||||||
|
2685.155556
|
||||||
|
2053.760000
|
||||||
|
3723.946667
|
||||||
|
3195.618462
|
||||||
|
2328.221538
|
||||||
|
1786.533333
|
||||||
|
7288.960000
|
||||||
|
4646.700000
|
||||||
|
2633.120000];
|
||||||
|
|
||||||
|
python = [2390.741
|
||||||
|
2597.527
|
||||||
|
2841.191
|
||||||
|
3066.652
|
||||||
|
3092.187
|
||||||
|
3444.710
|
||||||
|
3633.320
|
||||||
|
4333.845
|
||||||
|
5316.897
|
||||||
|
5201.486
|
||||||
|
4608.231
|
||||||
|
4481.357
|
||||||
|
3876.925
|
||||||
|
2961.753
|
||||||
|
2435.427
|
||||||
|
1344.871
|
||||||
|
606.953
|
||||||
|
298.559
|
||||||
|
120.772
|
||||||
|
50.369
|
||||||
|
17.033];
|
||||||
|
|
||||||
|
plot(log2(size), mflops,log2(size), dspl, log2(size), python)
|
|
@ -0,0 +1,31 @@
|
||||||
|
# -*- coding: utf-8 -*-
|
||||||
|
"""
|
||||||
|
Created on Tue Sep 22 14:48:55 2020
|
||||||
|
|
||||||
|
@author: bakhu
|
||||||
|
"""
|
||||||
|
|
||||||
|
import numpy as np
|
||||||
|
import time
|
||||||
|
import scipy.fftpack as fft
|
||||||
|
|
||||||
|
N = 4194304
|
||||||
|
|
||||||
|
x = np.random.randn(N) + 1j * np.random.randn(N)
|
||||||
|
|
||||||
|
|
||||||
|
k = 4
|
||||||
|
s = N
|
||||||
|
mflops = np.zeros(22)
|
||||||
|
for q in range(22):
|
||||||
|
z = np.copy(x[0:s])
|
||||||
|
t = time.time()
|
||||||
|
for i in range(k):
|
||||||
|
y = fft.fft(z)
|
||||||
|
|
||||||
|
dt = (time.time() - t)/float(k)
|
||||||
|
mflops[21-q] = 5.0 * s * np.log2(s) / dt / 1E6
|
||||||
|
|
||||||
|
print("size = %8d Mflops = %12.3f" % (s, mflops[21-q]))
|
||||||
|
k = int(k * 1.7)
|
||||||
|
s = int(s / 2)
|
|
@ -0,0 +1,113 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <time.h>
|
||||||
|
#include "dspl.h"
|
||||||
|
|
||||||
|
#define NMAX 4194304
|
||||||
|
#define L 18
|
||||||
|
#define SIZE_FACTOR 2.3
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
int fft_perf_cmplx(int* len, double* dlen, double* mflops)
|
||||||
|
{
|
||||||
|
int k = (int)( 8 * pow(SIZE_FACTOR, L));
|
||||||
|
int i, j, err;
|
||||||
|
clock_t t;
|
||||||
|
double dt;
|
||||||
|
complex_t *x = NULL;
|
||||||
|
complex_t *X = NULL;
|
||||||
|
|
||||||
|
x = (complex_t*) malloc (NMAX * sizeof(complex_t));
|
||||||
|
if(!x)
|
||||||
|
{
|
||||||
|
err = ERROR_MALLOC;
|
||||||
|
goto exit_label;
|
||||||
|
}
|
||||||
|
X = (complex_t*) malloc (NMAX * sizeof(complex_t));
|
||||||
|
if(!X)
|
||||||
|
{
|
||||||
|
err = ERROR_MALLOC;
|
||||||
|
goto exit_label;
|
||||||
|
}
|
||||||
|
|
||||||
|
fft_t pfft = {0}; /* fft structure */
|
||||||
|
random_t rnd = {0}; /* random structure */
|
||||||
|
|
||||||
|
err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
|
||||||
|
|
||||||
|
err = randn((double*)x, 2*NMAX, 0.0, 1.0, &rnd);
|
||||||
|
if(err != RES_OK)
|
||||||
|
goto exit_label;
|
||||||
|
|
||||||
|
|
||||||
|
printf("--------------------\n");
|
||||||
|
printf("FFT size MFlops\n");
|
||||||
|
printf("--------------------\n");
|
||||||
|
for(j = 0; j < L; j++)
|
||||||
|
{
|
||||||
|
dlen[j] = len[j];
|
||||||
|
t = clock();
|
||||||
|
/* FFT WORK CYCLE*/
|
||||||
|
for(i = 0; i < k; i++)
|
||||||
|
{
|
||||||
|
fft_cmplx(x, len[j], &pfft, X);
|
||||||
|
}
|
||||||
|
|
||||||
|
t = clock() - t;
|
||||||
|
|
||||||
|
dt = ((double)t)/CLOCKS_PER_SEC/(double)k;
|
||||||
|
|
||||||
|
mflops[j] = 5.0 * (double)len[j] *
|
||||||
|
log2((double)len[j]) / dt / 1E6;
|
||||||
|
|
||||||
|
printf ("%8d %10.1f (k = %d)\n", len[j], mflops[j], k);
|
||||||
|
|
||||||
|
k /= SIZE_FACTOR;
|
||||||
|
}
|
||||||
|
|
||||||
|
exit_label:
|
||||||
|
fft_free(&pfft);
|
||||||
|
free(x);
|
||||||
|
free(X);
|
||||||
|
return err;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
int main(int argc, char* argv[])
|
||||||
|
{
|
||||||
|
void* hdspl; /* DSPL handle */
|
||||||
|
void* hplot; /* GNUPLOT handles */
|
||||||
|
|
||||||
|
|
||||||
|
hdspl = dspl_load(); /* Load DSPL function */
|
||||||
|
|
||||||
|
int len_r2[L] = {2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096,
|
||||||
|
8192, 16384, 32768, 65536, 131072, 262144};
|
||||||
|
|
||||||
|
int len_nr[L] = {6, 9, 12, 15, 18, 24, 36, 80, 108, 210, 504, 1000,
|
||||||
|
1960, 4725, 10368, 27000, 75600, 165375};
|
||||||
|
|
||||||
|
int err;
|
||||||
|
double mflops[L] = {0};
|
||||||
|
double dlen[L];
|
||||||
|
|
||||||
|
printf("\n\nDouble precision complex 1D radix-2\n");
|
||||||
|
fft_perf_cmplx(len_r2, dlen, mflops);
|
||||||
|
writetxt(dlen, mflops, L, "dat/fft_cmplx_dspl_r2.txt");
|
||||||
|
|
||||||
|
printf("\n\nDouble precision complex 1D non-powers of two\n");
|
||||||
|
fft_perf_cmplx(len_nr, dlen, mflops);
|
||||||
|
writetxt(dlen, mflops, L, "dat/fft_cmplx_dspl_nr.txt");
|
||||||
|
|
||||||
|
|
||||||
|
exit_label:
|
||||||
|
/* free dspl handle */
|
||||||
|
dspl_free(hdspl);
|
||||||
|
return err;
|
||||||
|
}
|
||||||
|
|
|
@ -39,4 +39,5 @@ bin/$(LIB_NAME):$(RELEASE_DIR)/$(LIB_NAME)
|
||||||
clean:
|
clean:
|
||||||
rm -f $(OBJ_DIR)/*.o
|
rm -f $(OBJ_DIR)/*.o
|
||||||
rm -f bin/*.exe
|
rm -f bin/*.exe
|
||||||
|
rm -f bin/*.dll
|
||||||
|
rm -f bin/*.so
|
||||||
|
|
|
@ -0,0 +1,21 @@
|
||||||
|
clear all; close all; clc;
|
||||||
|
addpath('octave');
|
||||||
|
|
||||||
|
x = readbin('dat/x_fft_double_radix2.dat');
|
||||||
|
y = fft(x);
|
||||||
|
writebin(y, 1, 'dat/y_fft_double_radix2.dat');
|
||||||
|
|
||||||
|
|
||||||
|
x = readbin('dat/x_fft_complex_radix2.dat');
|
||||||
|
y = fft(x);
|
||||||
|
writebin(y, 1, 'dat/y_fft_complex_radix2.dat');
|
||||||
|
|
||||||
|
|
||||||
|
x = readbin('dat/x_fft_double_common.dat');
|
||||||
|
y = fft(x);
|
||||||
|
writebin(y, 1, 'dat/y_fft_double_common.dat');
|
||||||
|
|
||||||
|
|
||||||
|
x = readbin('dat/x_fft_complex_common.dat');
|
||||||
|
y = fft(x);
|
||||||
|
writebin(y, 1, 'dat/y_fft_complex_common.dat');
|
|
@ -0,0 +1,15 @@
|
||||||
|
clear all; close all; clc;
|
||||||
|
addpath('octave');
|
||||||
|
|
||||||
|
x = readbin('dat/real.dat');
|
||||||
|
y = mean(x);
|
||||||
|
writebin(y, 0, 'dat/mean_real.dat');
|
||||||
|
y = std(x);
|
||||||
|
writebin(y, 0, 'dat/std_real.dat');
|
||||||
|
|
||||||
|
x = readbin('dat/complex.dat');
|
||||||
|
y = mean(x);
|
||||||
|
writebin(y, 1, 'dat/mean_cmplx.dat');
|
||||||
|
y = std(x);
|
||||||
|
writebin(y, 0, 'dat/std_cmplx.dat');
|
||||||
|
|
|
@ -1,9 +1,10 @@
|
||||||
clear all; close all; clc;
|
clear all; close all; clc;
|
||||||
addpath('octave');
|
addpath('octave');
|
||||||
|
|
||||||
x = readbin('dat/x.dat');
|
x = readbin('dat/real.dat');
|
||||||
if(isreal(x))
|
writebin(x, 0, 'dat/yreal.dat');
|
||||||
writebin(x, 0, 'dat/y.dat');
|
|
||||||
else
|
x = readbin('dat/complex.dat');
|
||||||
writebin(x, 1, 'dat/y.dat');
|
writebin(x, 1, 'dat/ycomplex.dat');
|
||||||
end
|
|
||||||
|
|
||||||
|
|
|
@ -1,9 +0,0 @@
|
||||||
#ifndef DSPL_VERIF_H
|
|
||||||
#define DSPL_VERIF_H
|
|
||||||
|
|
||||||
#define VERIF_STR_LEN 48
|
|
||||||
#define VERIF_CHAR_POINT 46
|
|
||||||
#define VERIF_LEVEL_COMPLEX 1E-11
|
|
||||||
#define VERIF_LEVEL_DOUBLE 1E-12
|
|
||||||
|
|
||||||
#endif
|
|
|
@ -0,0 +1,76 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <time.h>
|
||||||
|
#include "dspl.h"
|
||||||
|
|
||||||
|
|
||||||
|
/* 2^15 */
|
||||||
|
#define FFT_R2_SIZE 32768
|
||||||
|
/* 2 * 3 * 5 * 7 * 11 * 13 */
|
||||||
|
#define FFT_CT_SIZE 30030
|
||||||
|
|
||||||
|
int main(int argc, char* argv[])
|
||||||
|
{
|
||||||
|
void* hdspl; /* DSPL handle */
|
||||||
|
fft_t pfft = {0};
|
||||||
|
|
||||||
|
int verr, nx, type;
|
||||||
|
double derr;
|
||||||
|
complex_t *yout = NULL;
|
||||||
|
complex_t *xc = NULL;
|
||||||
|
double *xd = NULL;
|
||||||
|
|
||||||
|
hdspl = dspl_load(); /* Load DSPL function */
|
||||||
|
|
||||||
|
verif_data_gen(FFT_R2_SIZE, DAT_DOUBLE, "dat/x_fft_double_radix2.dat");
|
||||||
|
verif_data_gen(FFT_R2_SIZE, DAT_COMPLEX, "dat/x_fft_complex_radix2.dat");
|
||||||
|
verif_data_gen(FFT_CT_SIZE, DAT_DOUBLE, "dat/x_fft_double_common.dat");
|
||||||
|
verif_data_gen(FFT_CT_SIZE, DAT_COMPLEX, "dat/x_fft_complex_common.dat");
|
||||||
|
|
||||||
|
yout = (complex_t*)malloc(FFT_R2_SIZE * sizeof(complex_t ));
|
||||||
|
|
||||||
|
system("octave octave/fft_verification.m");
|
||||||
|
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
readbin("dat/x_fft_double_radix2.dat", (void**)(&xd), &nx, &type);
|
||||||
|
fft(xd, FFT_R2_SIZE, &pfft, yout);
|
||||||
|
verif_str_cmplx(yout, FFT_R2_SIZE, "fft radix-2 for double dat",
|
||||||
|
"dat/y_fft_double_radix2.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
readbin("dat/x_fft_complex_radix2.dat", (void**)(&xc), &nx, &type);
|
||||||
|
fft_cmplx(xc, FFT_R2_SIZE, &pfft, yout);
|
||||||
|
verif_str_cmplx(yout, FFT_R2_SIZE, "fft radix-2 for complex dat",
|
||||||
|
"dat/y_fft_complex_radix2.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
readbin("dat/x_fft_double_common.dat", (void**)(&xd), &nx, &type);
|
||||||
|
fft(xd, FFT_CT_SIZE, &pfft, yout);
|
||||||
|
verif_str_cmplx(yout, FFT_CT_SIZE, "fft composite for double dat",
|
||||||
|
"dat/y_fft_double_common.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
readbin("dat/x_fft_complex_common.dat", (void**)(&xc), &nx, &type);
|
||||||
|
fft_cmplx(xc, FFT_CT_SIZE, &pfft, yout);
|
||||||
|
verif_str_cmplx(yout, FFT_CT_SIZE, "fft composite for complex dat",
|
||||||
|
"dat/y_fft_complex_common.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
/* free dspl handle */
|
||||||
|
dspl_free(hdspl);
|
||||||
|
|
||||||
|
if(yout)
|
||||||
|
free(yout);
|
||||||
|
if(xc)
|
||||||
|
free(xc);
|
||||||
|
if(xd)
|
||||||
|
free(xd);
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
|
@ -0,0 +1,66 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <time.h>
|
||||||
|
#include "dspl.h"
|
||||||
|
|
||||||
|
|
||||||
|
#define SIZE 500
|
||||||
|
|
||||||
|
|
||||||
|
int main(int argc, char* argv[])
|
||||||
|
{
|
||||||
|
void* hdspl; /* DSPL handle */
|
||||||
|
|
||||||
|
int err, verr, nx, type;
|
||||||
|
double derr;
|
||||||
|
complex_t yc;
|
||||||
|
double yd;
|
||||||
|
complex_t *xc = NULL;
|
||||||
|
double *xd = NULL;
|
||||||
|
|
||||||
|
|
||||||
|
hdspl = dspl_load(); /* Load DSPL function */
|
||||||
|
|
||||||
|
verif_data_gen(SIZE, DAT_DOUBLE, "dat/real.dat");
|
||||||
|
verif_data_gen(SIZE, DAT_COMPLEX, "dat/complex.dat");
|
||||||
|
|
||||||
|
system("octave octave/statistic_verification.m");
|
||||||
|
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
readbin("dat/real.dat", (void**)(&xd), &nx, &type);
|
||||||
|
mean(xd, SIZE, &yd);
|
||||||
|
verif_str(&yd, 1, "mean for double data:",
|
||||||
|
"dat/mean_real.dat",
|
||||||
|
"verification.log");
|
||||||
|
std(xd, SIZE, &yd);
|
||||||
|
verif_str(&yd, 1, "std for double data:", "dat/std_real.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
readbin("dat/complex.dat", (void**)(&xc), &nx, &type);
|
||||||
|
mean_cmplx(xc, SIZE, &yc);
|
||||||
|
verif_str_cmplx(&yc, 1, "mean for complex data:",
|
||||||
|
"dat/mean_cmplx.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
std_cmplx(xc, SIZE, &yd);
|
||||||
|
verif_str(&yd, 1, "std for complex data:",
|
||||||
|
"dat/std_cmplx.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
|
||||||
|
/* free dspl handle */
|
||||||
|
dspl_free(hdspl);
|
||||||
|
if(xc)
|
||||||
|
free(xc);
|
||||||
|
if(xd)
|
||||||
|
free(xd);
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
|
@ -0,0 +1,66 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <time.h>
|
||||||
|
#include "dspl.h"
|
||||||
|
|
||||||
|
#define SIZE 10000
|
||||||
|
|
||||||
|
int main(int argc, char* argv[])
|
||||||
|
{
|
||||||
|
void* hdspl; /* DSPL handle */
|
||||||
|
complex_t *xc = NULL;
|
||||||
|
double *xd = NULL;
|
||||||
|
|
||||||
|
int nx, type, err, verr;
|
||||||
|
double derr;
|
||||||
|
|
||||||
|
|
||||||
|
hdspl = dspl_load(); /* Load DSPL function */
|
||||||
|
|
||||||
|
/* genreation real random data for verification */
|
||||||
|
verif_data_gen(SIZE, DAT_DOUBLE, "dat/real.dat");
|
||||||
|
|
||||||
|
/* genreation complex random data for verification */
|
||||||
|
verif_data_gen(SIZE, DAT_COMPLEX, "dat/complex.dat");
|
||||||
|
|
||||||
|
/* RUN verification in octave */
|
||||||
|
system("octave octave/writebin_readbin_verification.m");
|
||||||
|
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
/* Read real input data from the file */
|
||||||
|
readbin("dat/real.dat", (void**)(&xd), &nx, &type);
|
||||||
|
|
||||||
|
/* Processing */
|
||||||
|
|
||||||
|
/* verification libdspl output and octave output */
|
||||||
|
verif_str(xd, SIZE, "writebin and readbin for real data:",
|
||||||
|
"dat/yreal.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------*/
|
||||||
|
/* Read complex input data from the file */
|
||||||
|
readbin("dat/complex.dat", (void**)(&xc), &nx, &type);
|
||||||
|
|
||||||
|
/* Processing */
|
||||||
|
|
||||||
|
/* verification libdspl output and octave output */
|
||||||
|
verif_str_cmplx(xc, SIZE, "writebin and readbin for complex data:",
|
||||||
|
"dat/ycomplex.dat",
|
||||||
|
"verification.log");
|
||||||
|
|
||||||
|
|
||||||
|
/* free dspl handle */
|
||||||
|
dspl_free(hdspl);
|
||||||
|
|
||||||
|
|
||||||
|
if(xc)
|
||||||
|
free(xc);
|
||||||
|
if(xd)
|
||||||
|
free(xd);
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
|
@ -1,96 +0,0 @@
|
||||||
#include <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <time.h>
|
|
||||||
#include "dspl.h"
|
|
||||||
#include "dspl_verif.h"
|
|
||||||
|
|
||||||
#define ARRAY_MAX_SIZE 1000000
|
|
||||||
|
|
||||||
|
|
||||||
int main(int argc, char* argv[])
|
|
||||||
{
|
|
||||||
void* hdspl; /* DSPL handle */
|
|
||||||
complex_t *pxd = NULL;
|
|
||||||
complex_t *pyd = NULL;
|
|
||||||
char str[512] = {0};
|
|
||||||
|
|
||||||
int nx, ny, tx, err, verr;
|
|
||||||
double derr;
|
|
||||||
|
|
||||||
random_t rnd = {0}; /* random structure */
|
|
||||||
|
|
||||||
hdspl = dspl_load(); /* Load DSPL function */
|
|
||||||
|
|
||||||
sprintf(str, "writebin and readbin for complex data:");
|
|
||||||
while(strlen(str) < VERIF_STR_LEN)
|
|
||||||
str[strlen(str)] = VERIF_CHAR_POINT;
|
|
||||||
|
|
||||||
|
|
||||||
err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
|
|
||||||
if(err != RES_OK)
|
|
||||||
goto exit_error_code;
|
|
||||||
|
|
||||||
err = randi(&nx, 1, 1, ARRAY_MAX_SIZE, &rnd);
|
|
||||||
if(err != RES_OK)
|
|
||||||
goto exit_error_code;
|
|
||||||
|
|
||||||
pxd = (complex_t*) malloc(nx * sizeof(complex_t));
|
|
||||||
|
|
||||||
|
|
||||||
/****************** verifiсation function start **************************/
|
|
||||||
err = randn((double*)pxd, 2*nx, 1.0, 10.0, &rnd);
|
|
||||||
if(err != RES_OK)
|
|
||||||
goto exit_error_code;
|
|
||||||
|
|
||||||
/************ Write input verification data to the file ******************/
|
|
||||||
err = writebin((void*)pxd, nx, DAT_COMPLEX, "dat/x.dat");
|
|
||||||
if(err != RES_OK)
|
|
||||||
goto exit_error_code;
|
|
||||||
|
|
||||||
/************ RUN external verificator (octave or python) ****************/
|
|
||||||
err = system("octave octave/writebin_readbin_verification.m");
|
|
||||||
|
|
||||||
/***************** Read external verificator output **********************/
|
|
||||||
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;
|
|
||||||
}
|
|
||||||
|
|
||||||
/**************************** Verification *******************************/
|
|
||||||
verr = verif_cmplx(pxd, pyd, ny, VERIF_LEVEL_COMPLEX, &derr);
|
|
||||||
if(verr != DSPL_VERIF_SUCCESS)
|
|
||||||
goto exit_error_verif;
|
|
||||||
|
|
||||||
sprintf(str, "%s ok (err = %12.4E)", str, derr);
|
|
||||||
goto exit_label;
|
|
||||||
|
|
||||||
exit_error_code:
|
|
||||||
sprintf(str, "%s FAILED (with code = 0x%8x)", str, err);
|
|
||||||
goto exit_label;
|
|
||||||
|
|
||||||
exit_error_verif:
|
|
||||||
sprintf(str, "%s FAILED (err = %12.4E)", str, derr);
|
|
||||||
|
|
||||||
exit_label:
|
|
||||||
/************************ write str to log file **************************/
|
|
||||||
|
|
||||||
printf("%s\n", str);
|
|
||||||
addlog(str, "verification.log");
|
|
||||||
|
|
||||||
if(pxd)
|
|
||||||
free(pxd);
|
|
||||||
if(pyd)
|
|
||||||
free(pyd);
|
|
||||||
|
|
||||||
/* free dspl handle */
|
|
||||||
dspl_free(hdspl);
|
|
||||||
|
|
||||||
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
|
@ -1,83 +0,0 @@
|
||||||
#include <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <time.h>
|
|
||||||
#include "dspl.h"
|
|
||||||
#include "dspl_verif.h"
|
|
||||||
|
|
||||||
#define ARRAY_MAX_SIZE 1000000
|
|
||||||
|
|
||||||
int main(int argc, char* argv[])
|
|
||||||
{
|
|
||||||
void* hdspl; /* DSPL handle */
|
|
||||||
double* pxd = NULL;
|
|
||||||
double* pyd = NULL;
|
|
||||||
char str[512] = {0};
|
|
||||||
|
|
||||||
|
|
||||||
int nx, ny, tx, err, verr;
|
|
||||||
double derr;
|
|
||||||
|
|
||||||
random_t rnd = {0}; /* random structure */
|
|
||||||
|
|
||||||
hdspl = dspl_load(); /* Load DSPL function */
|
|
||||||
|
|
||||||
sprintf(str, "writebin and readbin for double data:");
|
|
||||||
while(strlen(str) < VERIF_STR_LEN)
|
|
||||||
str[strlen(str)] = VERIF_CHAR_POINT;
|
|
||||||
|
|
||||||
|
|
||||||
err = random_init(&rnd, RAND_TYPE_MRG32K3A, NULL);
|
|
||||||
if(err != RES_OK)
|
|
||||||
goto exit_error_code;
|
|
||||||
|
|
||||||
err = randi(&nx, 1, 1, ARRAY_MAX_SIZE, &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, VERIF_LEVEL_DOUBLE, &derr);
|
|
||||||
if(verr != DSPL_VERIF_SUCCESS)
|
|
||||||
goto exit_error_verif;
|
|
||||||
|
|
||||||
sprintf(str, "%s ok (err = %12.4E)", str, derr);
|
|
||||||
goto exit_label;
|
|
||||||
|
|
||||||
exit_error_code:
|
|
||||||
sprintf(str, "%s FAILED (with code = 0x%8x)", str, err);
|
|
||||||
goto exit_label;
|
|
||||||
exit_error_verif:
|
|
||||||
sprintf(str, "%s FAILED (err = %12.4E)", str, derr);
|
|
||||||
|
|
||||||
exit_label:
|
|
||||||
addlog(str, "verification.log");
|
|
||||||
printf("%s\n", str);
|
|
||||||
if(pxd)
|
|
||||||
free(pxd);
|
|
||||||
if(pyd)
|
|
||||||
free(pyd);
|
|
||||||
/* free dspl handle */
|
|
||||||
dspl_free(hdspl);
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
Ładowanie…
Reference in New Issue