added doxydoc for fft_cmplx and ifft_cmplx_test

added functions histogram, histogram_norm and minmax
pull/6/merge
Dsplib 2019-04-03 23:37:47 +03:00
rodzic 163dd4d4e8
commit 90a32fdd21
12 zmienionych plików z 432 dodań i 57 usunięć

Wyświetl plik

@ -2,6 +2,7 @@
body, table, div, p, dl {
font: 400 14px/22px Roboto,sans-serif;
background-color: #EFEFEF;
}
p.reference, p.definition {

Wyświetl plik

@ -19,18 +19,17 @@ $generatedby &#160;<a href="http://www.doxygen.org/index.html">
<!--END !GENERATE_TREEVIEW-->
</div> <!--div id="dsplib_content"-->
</div><!--div id="dsplib_maket"-->
</div> <!-- div id="dsplib_main" -->
</div> <!-- div id="wrapper" -->
<hr align="center" width="100%" size="1" color="#000000" />
<div class="dsplib_copyright">
&copy; Бахурин Сергей 2015 - 2019. Все права защищены.
Любое копирование материалов сайте без разрешения автора запрещено.
</div>
</div> <!-- div id="wrapper" -->
</body>
</html>

Wyświetl plik

@ -98,7 +98,7 @@ $extrastylesheet
<div id="dsplib_content">

Wyświetl plik

@ -91,6 +91,86 @@ www.dsplib.org
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn int ifft_cmplx(complex_t* x, int n, fft_t* pfft, complex_t* y)
\brief Обратное быстрое преобразование Фурье
Функция рассчитывает \f$ n \f$-точечное обратное быстрое преобразование Фурье
от \f$ x(m) \f$, \f$ m = 0 \ldots n-1 \f$.<BR>
\f[
Y(k) = \frac{1}{N} \sum_{m = 0}^{n-1} x(m) \exp
\left( j \frac{2\pi}{n} m k \right),
\f]
где \f$ k = 0 \ldots n-1 \f$.
Для расчета используется алгоритм БПФ составной длины.
\param[in] x Указатель на входной комплексный вектор \f$x(m)\f$,
\f$ m = 0 \ldots n-1 \f$. <BR>
Размер вектора `[n x 1]`. <BR><BR>
\param[in] n Размер ОБПФ \f$n\f$.<BR>
Размер ОБПФ может быть составным вида
\f$n = n_0 \times n_1 \times n_2 \times n_3 \times \ldots \times n_p \times m\f$,
где \f$n_i = 2,3,5,7\f$, а \f$m \f$ --
произвольный простой множитель не превосходящий 46340
(см. описание функции \ref fft_create).
<BR><BR>
\param[in] pfft Указатель на структуру `fft_t`. <BR>
Указатель не должен быть `NULL`. <BR>
Структура \ref fft_t должна быть предварительно однократно
заполнена функцией \ref fft_create, и память должна быть
очищена перед выходом функцией \ref fft_free.
<BR><BR>
\param[out] y Указатель на вектор результата ОБПФ \f$Y(k)\f$,
\f$ k = 0 \ldots n-1 \f$.
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` если расчет произведен успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".<BR><BR>
Пример использования функции `fft`:
\include ifft_cmplx_test.c
Результат работы программы:
\verbatim
| x[ 0] = 1.000 0.000 | y[ 0] = -0.517 0.686 | z[ 0] = 1.000 0.000 |
| x[ 1] = 0.540 0.841 | y[ 1] = -0.943 0.879 | z[ 1] = 0.540 0.841 |
| x[ 2] = -0.416 0.909 | y[ 2] = -2.299 1.492 | z[ 2] = -0.416 0.909 |
| x[ 3] = -0.990 0.141 | y[ 3] = 16.078 -6.820 | z[ 3] = -0.990 0.141 |
| x[ 4] = -0.654 -0.757 | y[ 4] = 2.040 -0.470 | z[ 4] = -0.654 -0.757 |
| x[ 5] = 0.284 -0.959 | y[ 5] = 1.130 -0.059 | z[ 5] = 0.284 -0.959 |
| x[ 6] = 0.960 -0.279 | y[ 6] = 0.786 0.097 | z[ 6] = 0.960 -0.279 |
| x[ 7] = 0.754 0.657 | y[ 7] = 0.596 0.183 | z[ 7] = 0.754 0.657 |
| x[ 8] = -0.146 0.989 | y[ 8] = 0.470 0.240 | z[ 8] = -0.146 0.989 |
| x[ 9] = -0.911 0.412 | y[ 9] = 0.375 0.283 | z[ 9] = -0.911 0.412 |
| x[10] = -0.839 -0.544 | y[10] = 0.297 0.318 | z[10] = -0.839 -0.544 |
| x[11] = 0.004 -1.000 | y[11] = 0.227 0.350 | z[11] = 0.004 -1.000 |
| x[12] = 0.844 -0.537 | y[12] = 0.161 0.380 | z[12] = 0.844 -0.537 |
| x[13] = 0.907 0.420 | y[13] = 0.094 0.410 | z[13] = 0.907 0.420 |
| x[14] = 0.137 0.991 | y[14] = 0.023 0.442 | z[14] = 0.137 0.991 |
| x[15] = -0.760 0.650 | y[15] = -0.059 0.479 | z[15] = -0.760 0.650 |
| x[16] = -0.958 -0.288 | y[16] = -0.161 0.525 | z[16] = -0.958 -0.288 |
| x[17] = -0.275 -0.961 | y[17] = -0.300 0.588 | z[17] = -0.275 -0.961 |
\endverbatim
\author
Бахурин Сергей.
www.dsplib.org
***************************************************************************** */
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn int fft(double* x, int n, fft_t* pfft, complex_t* y)
@ -166,9 +246,92 @@ www.dsplib.org
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn int fft_create(fft_t *pfft, int n)
\fn int fft_cmplx(complex_t* x, int n, fft_t* pfft, complex_t* y)
\brief Быстрое преобразование Фурье комплексного сигнала
Функция рассчитывает \f$ n \f$-точечное быстрое преобразование Фурье
комплексного сигнала \f$ x(m) \f$, \f$ m = 0 \ldots n-1 \f$.<BR>
\f[
Y(k) = \sum_{m = 0}^{n-1} x(m) \exp
\left( -j \frac{2\pi}{n} m k \right),
\f]
где \f$ k = 0 \ldots n-1 \f$.
Для расчета используется алгоритм БПФ составной длины.
\param[in] x Указатель на вектор комплексного
входного сигнала \f$x(m)\f$,
\f$ m = 0 \ldots n-1 \f$. <BR>
Размер вектора `[n x 1]`. <BR><BR>
\param[in] n Размер БПФ \f$n\f$.<BR>
Размер БПФ может быть составным вида
\f$n = n_0 \times n_1 \times n_2 \times n_3 \times \ldots \times n_p \times m\f$,
где \f$n_i = 2,3,5,7\f$, а \f$m \f$ --
произвольный простой множитель не превосходящий 46340
(см. описание функции \ref fft_create).
<BR><BR>
\param[in] pfft Указатель на структуру `fft_t`. <BR>
Указатель не должен быть `NULL`. <BR>
Структура \ref fft_t должна быть предварительно однократно
заполнена функцией \ref fft_create, и память должна быть
очищена перед выходом функцией \ref fft_free.
<BR><BR>
\param[out] y Указатель на комплексный вектор
результата БПФ \f$Y(k)\f$,
\f$ k = 0 \ldots n-1 \f$.
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` если расчет произведен успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".<BR><BR>
Пример использования функции `fft`:
\include fft_cmplx_test.c
Результат работы программы:
\verbatim
y[ 0] = -0.517 0.686
y[ 1] = -0.943 0.879
y[ 2] = -2.299 1.492
y[ 3] = 16.078 -6.820
y[ 4] = 2.040 -0.470
y[ 5] = 1.130 -0.059
y[ 6] = 0.786 0.097
y[ 7] = 0.596 0.183
y[ 8] = 0.470 0.240
y[ 9] = 0.375 0.283
y[10] = 0.297 0.318
y[11] = 0.227 0.350
y[12] = 0.161 0.380
y[13] = 0.094 0.410
y[14] = 0.023 0.442
y[15] = -0.059 0.479
y[16] = -0.161 0.525
y[17] = -0.300 0.588
\endverbatim
\author
Бахурин Сергей.
www.dsplib.org
***************************************************************************** */
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn int fft_create(fft_t* pfft, int n)
\brief Заполнение структуры `fft_t` для алгоритма БПФ
Функция производит выделение памяти и рассчет векторов
@ -179,29 +342,29 @@ www.dsplib.org
\param[in,out] pfft Указатель на структуру `fft_t`. <BR>
Указатель не должен быть `NULL`. <BR><BR>
\param[in] n Размер БПФ \f$n\f$.<BR>
Размер БПФ может быть составным вида
\f$n = n_0 \times n_1 \times n_2 \times n_3 \times \ldots \times n_p \times m\f$,
где \f$n_i = 2,3,5,7\f$, а \f$m \f$ --
произвольный простой множитель не превосходящий 46340.
<BR><BR>
Таким образом алгоритм БПФ поддерживает произвольные
длины, равные целой степени чисел 2,3,5,7,
а также различные их комбинации. <BR><BR>
Так например, при \f$ n = 725760 \f$ структура
будет успешно заполнена, отому что
\f$725760 = 2 \cdot 3 \cdot 4 \cdot 5 \cdot 6 \cdot 7 \cdot 9 \cdot 16 \f$,
т.е. получается как произведение множителей 2,3,5,7.
<BR><BR>
При \f$ n = 172804 = 43201 \cdot 4 \f$ структура также
будет успешно заполнена, потому что простой множитель
входящий в \f$n\f$ не превосходит 46340.
<BR><BR>
Для размера
\f$ n = 13 \cdot 17 \cdot 23 \cdot 13 = 66079 \f$
функция вернет ошибку, поскольку 66079 больше 46340
и не является результатом произведения чисел 2,3,5,7.
<BR><BR>
\param[in] n Размер БПФ \f$n\f$.<BR>
Размер БПФ может быть составным вида
\f$n = n_0 \times n_1 \times n_2 \times n_3 \times \ldots \times n_p \times m\f$,
где \f$n_i = 2,3,5,7\f$, а \f$m \f$ --
произвольный простой множитель не превосходящий 46340.
<BR><BR>
Таким образом алгоритм БПФ поддерживает произвольные
длины, равные целой степени чисел 2,3,5,7,
а также различные их комбинации. <BR><BR>
Так например, при \f$ n = 725760 \f$ структура
будет успешно заполнена, отому что
\f$725760 = 2 \cdot 3 \cdot 4 \cdot 5 \cdot 6 \cdot 7 \cdot 9 \cdot 16 \f$,
т.е. получается как произведение множителей 2,3,5,7.
<BR><BR>
При \f$ n = 172804 = 43201 \cdot 4 \f$ структура также
будет успешно заполнена, потому что простой множитель
входящий в \f$n\f$ не превосходит 46340.
<BR><BR>
Для размера
\f$ n = 13 \cdot 17 \cdot 23 \cdot 13 = 66079 \f$
функция вернет ошибку, поскольку 66079 больше 46340
и не является результатом произведения чисел 2,3,5,7.
<BR><BR>
\return
`RES_OK` если структура заполнена успешно. <BR>

Wyświetl plik

@ -64,7 +64,7 @@ int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
/*******************************************************************************
Real vector FFT
*******************************************************************************/
int DSPL_API fft(double *x, int n, fft_t* pfft, complex_t* y)
int DSPL_API fft(double* x, int n, fft_t* pfft, complex_t* y)
{
int err;
@ -89,7 +89,7 @@ int DSPL_API fft(double *x, int n, fft_t* pfft, complex_t* y)
/*******************************************************************************
COMPLEX vector FFT
*******************************************************************************/
int DSPL_API fft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
int DSPL_API fft_cmplx(complex_t* x, int n, fft_t* pfft, complex_t* y)
{
int err;
@ -205,7 +205,7 @@ label_size:
/*******************************************************************************
FFT create for composite N
*******************************************************************************/
int DSPL_API fft_create(fft_t *pfft, int n)
int DSPL_API fft_create(fft_t* pfft, int n)
{
int n1, n2, addr, s, k, m, nw, err;

Wyświetl plik

@ -0,0 +1,122 @@
/*
* 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 DSPL_API histogram(double* x, int n, int nh, double* pedges, double* ph)
{
double xmin, xmax;
int k, ind;
int res;
if(!x || !pedges || !ph)
return ERROR_PTR;
if(n<1 || nh<1)
return ERROR_SIZE;
res = minmax(x, n, &xmin, &xmax);
if(res != RES_OK)
return res;
res = linspace(xmin, xmax, nh+1, DSPL_SYMMETRIC, pedges);
if(res != RES_OK)
return res;
memset(ph, 0, nh*sizeof(double));
for(k = 0; k < n; k++)
{
ind = 0;
while(ind<nh && x[k]>=pedges[ind])
ind++;
ph[ind-1]+=1.0;
}
return RES_OK;
}
int DSPL_API histogram_norm(double* y, int n, int nh, double* x, double* w)
{
double *pedges = NULL;
int k, res;
if(!y || !x || !w)
return ERROR_PTR;
if(n<1 || nh<1)
return ERROR_SIZE;
pedges = (double*)malloc((nh+1)*sizeof(double));
res = histogram(y, n, nh, pedges, w);
if(res != RES_OK)
goto exit_label;
for(k = 1; k < nh+1; k++)
{
x[k-1] = 0.5*(pedges[k] + pedges[k-1]);
w[k-1] /= ((double)n * (pedges[k] - pedges[k-1]));
}
res = RES_OK;
exit_label:
if(pedges)
free(pedges);
return res;
}
int DSPL_API minmax(double* x, int n, double* xmin, double* xmax)
{
int k;
double min, max;
if(!x)
return ERROR_PTR;
if(n<1)
return ERROR_SIZE;
min = max = x[0];
for(k = 1; k < n; k++)
{
min = x[k] < min ? x[k] : min;
max = x[k] > max ? x[k] : max;
}
if(xmin)
*xmin = min;
if(xmax)
*xmax = max;
return RES_OK;
}

Wyświetl plik

@ -3,37 +3,37 @@
#include <string.h>
#include "dspl.h"
#define N 16
#define N 18
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
void* handle; // DSPL handle
handle = dspl_load(); // Загрузка DSPL
complex_t x[N]; // complex input signal
complex_t y[N]; // DFT
fft_t pfft; // FFT object
//clear fft object
memset(&pfft, 0, sizeof(fft_t));
// Create FFT object
fft_create(&pfft, N);
for(int k = 0; k < N; k++)
{
RE(x[k]) = (double)k;
IM(x[k]) = 0.0;
}
//FFT
fft_cmplx(x, N, &pfft, y);
complex_t x[N]; // массив входного сигнала
complex_t y[N]; // массив результата БПФ
fft_t pfft; // FFT объект
memset(&pfft, 0, sizeof(fft_t)); // Заполняем FFT структуру нулями
fft_create(&pfft, N); // Создаем FFT структуру для длины N
// заполняем массив входного сигнала
for(int k = 0; k < N; k++)
{
RE(x[k]) = (double)cos((double)k);
IM(x[k]) = (double)sin((double)k);
}
fft_cmplx(x, N, &pfft, y); // FFT
for(int k = 0; k < N; k++)
printf("y[%2d] = %9.3f%9.3f\n", k, RE(y[k]), IM(y[k]));
// Печать результата
for(int k = 0; k < N; k++)
printf("y[%2d] = %9.3f%9.3f\n", k, RE(y[k]), IM(y[k]));
fft_free(&pfft); // clear FFT object
dspl_free(handle); // free dspl handle
return 0;
fft_free(&pfft); // Очищаем структуру fft_t
dspl_free(handle); // Очищаем dspl handle
return 0;
}

Wyświetl plik

@ -0,0 +1,44 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 18
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Загрузка DSPL
complex_t x[N]; // массив входного сигнала
complex_t y[N]; // массив результата БПФ
complex_t z[N]; // массив результата ОБПФ
fft_t pfft; // FFT объект
memset(&pfft, 0, sizeof(fft_t)); // Заполняем FFT структуру нулями
fft_create(&pfft, N); // Создаем FFT структуру для длины N
// заполняем массив входного сигнала
for(int k = 0; k < N; k++)
{
RE(x[k]) = (double)cos((double)k);
IM(x[k]) = (double)sin((double)k);
}
fft_cmplx(x, N, &pfft, y); // FFT
ifft_cmplx(y, N, &pfft, z); // IFFT
// Печать результата
for(int k = 0; k < N; k++)
{
printf("| x[%2d] = %9.3f%9.3f ", k, RE(x[k]), IM(x[k]));
printf("| y[%2d] = %9.3f%9.3f ", k, RE(y[k]), IM(y[k]));
printf("| z[%2d] = %9.3f%9.3f |\n", k, RE(z[k]), IM(z[k]));
}
fft_free(&pfft); // Очищаем структуру fft_t
dspl_free(handle); // Очищаем dspl handle
return 0;
}

Wyświetl plik

@ -98,6 +98,8 @@ p_freqs2time freqs2time ;
p_freqz freqz ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_histogram histogram ;
p_histogram_norm histogram_norm ;
p_linspace linspace ;
p_log_cmplx log_cmplx ;
p_logspace logspace ;
@ -112,6 +114,7 @@ p_matrix_swap matrix_swap ;
p_matrix_swap_rows matrix_swap_rows ;
p_matrix_transpose matrix_transpose ;
p_matrix_transpose_hermite matrix_transpose_hermite ;
p_minmax minmax ;
p_poly_z2a_cmplx poly_z2a_cmplx ;
p_polyval polyval ;
p_polyval_cmplx polyval_cmplx ;
@ -254,6 +257,8 @@ void* dspl_load()
LOAD_FUNC(freqs2time);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(histogram);
LOAD_FUNC(histogram_norm);
LOAD_FUNC(linspace);
LOAD_FUNC(log_cmplx);
LOAD_FUNC(logspace);
@ -268,6 +273,7 @@ void* dspl_load()
LOAD_FUNC(matrix_swap_rows);
LOAD_FUNC(matrix_transpose);
LOAD_FUNC(matrix_transpose_hermite);
LOAD_FUNC(minmax);
LOAD_FUNC(poly_z2a_cmplx);
LOAD_FUNC(polyval);
LOAD_FUNC(polyval_cmplx);

Wyświetl plik

@ -588,6 +588,18 @@ DECLARE_FUNC(int, goertzel_cmplx, complex_t*
COMMA int
COMMA complex_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, histogram, double* x
COMMA int n
COMMA int nh
COMMA double* pedges
COMMA double* ph);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, histogram_norm, double* y
COMMA int n
COMMA int nh
COMMA double* x
COMMA double* w);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, linspace, double
COMMA double
COMMA int
@ -665,6 +677,11 @@ DECLARE_FUNC(int, matrix_transpose, matrix_t* a
DECLARE_FUNC(int, matrix_transpose_hermite, matrix_t* a
COMMA matrix_t* b);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, minmax, double* x
COMMA int n
COMMA double* xmin
COMMA double* xmax);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, poly_z2a_cmplx, complex_t*
COMMA int
COMMA int

Wyświetl plik

@ -98,6 +98,8 @@ p_freqs2time freqs2time ;
p_freqz freqz ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_histogram histogram ;
p_histogram_norm histogram_norm ;
p_linspace linspace ;
p_log_cmplx log_cmplx ;
p_logspace logspace ;
@ -112,6 +114,7 @@ p_matrix_swap matrix_swap ;
p_matrix_swap_rows matrix_swap_rows ;
p_matrix_transpose matrix_transpose ;
p_matrix_transpose_hermite matrix_transpose_hermite ;
p_minmax minmax ;
p_poly_z2a_cmplx poly_z2a_cmplx ;
p_polyval polyval ;
p_polyval_cmplx polyval_cmplx ;
@ -254,6 +257,8 @@ void* dspl_load()
LOAD_FUNC(freqs2time);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(histogram);
LOAD_FUNC(histogram_norm);
LOAD_FUNC(linspace);
LOAD_FUNC(log_cmplx);
LOAD_FUNC(logspace);
@ -268,6 +273,7 @@ void* dspl_load()
LOAD_FUNC(matrix_swap_rows);
LOAD_FUNC(matrix_transpose);
LOAD_FUNC(matrix_transpose_hermite);
LOAD_FUNC(minmax);
LOAD_FUNC(poly_z2a_cmplx);
LOAD_FUNC(polyval);
LOAD_FUNC(polyval_cmplx);

Wyświetl plik

@ -588,6 +588,18 @@ DECLARE_FUNC(int, goertzel_cmplx, complex_t*
COMMA int
COMMA complex_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, histogram, double* x
COMMA int n
COMMA int nh
COMMA double* pedges
COMMA double* ph);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, histogram_norm, double* y
COMMA int n
COMMA int nh
COMMA double* x
COMMA double* w);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, linspace, double
COMMA double
COMMA int
@ -665,6 +677,11 @@ DECLARE_FUNC(int, matrix_transpose, matrix_t* a
DECLARE_FUNC(int, matrix_transpose_hermite, matrix_t* a
COMMA matrix_t* b);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, minmax, double* x
COMMA int n
COMMA double* xmin
COMMA double* xmax);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, poly_z2a_cmplx, complex_t*
COMMA int
COMMA int