diff --git a/dox/doxy_stylesheet.css b/dox/doxy_stylesheet.css index 9c8d13a..5434f50 100644 --- a/dox/doxy_stylesheet.css +++ b/dox/doxy_stylesheet.css @@ -1,8 +1,15 @@ /* The standard CSS for doxygen 1.8.13 */ +/* body, table, div, p, dl { font: 400 14px/22px Roboto,sans-serif; } +*/ + +code{ + font-size:12pt; + font-weight: bold; +} p.reference, p.definition { font: 400 14px/22px Roboto,sans-serif; diff --git a/dox/footer.html b/dox/footer.html index b7b2b6f..1e17f0f 100644 --- a/dox/footer.html +++ b/dox/footer.html @@ -18,18 +18,7 @@ $generatedby   + - - - - - - - -
-© Бахурин Сергей 2015 - 2019. Все права защищены. -Любое копирование материалов сайте без разрешения автора запрещено. -
- - + \ No newline at end of file diff --git a/dox/header.html b/dox/header.html index 7c9ab13..e86cbf7 100644 --- a/dox/header.html +++ b/dox/header.html @@ -20,11 +20,16 @@ $treeview $search $mathjax - + + + + + + $extrastylesheet - + @@ -64,36 +69,27 @@ $extrastylesheet - -
-
- -
-
- +
+
+
- - - -
-
- - -
+ +
diff --git a/dspl/dox/ru/conv.dox b/dspl/dox/ru/conv.dox index 78283d7..47a9f9d 100644 --- a/dspl/dox/ru/conv.dox +++ b/dspl/dox/ru/conv.dox @@ -134,6 +134,100 @@ cc[5] = 0.0+24.0j +/*! **************************************************************************** +\ingroup FILTER_CONV_GROUP +\fn int conv_fft(double* a, int na, double* b, int nb, + fft_t* pfft, double* c) +\brief Линейная свертка двух вещественных векторов с использованием алгоритмов +быстрого преобразования Фурье + +Функция рассчитывает линейную свертку двух векторов \f$ c = a * b\f$ используя +секционную обработку с перекрытием в частотной области. Это позволяет сократить +вычислительные операции при расчете длинных сверток. + + + +\param[in] a Указатель на первый вектор \f$a\f$.
+ Размер вектора `[na x 1]`.

+ +\param[in] na Размер первого вектора.

+ +\param[in] b Указатель на второй вектор \f$b\f$.
+ Размер вектора `[nb x 1]`.

+ +\param[in] nb Размер второго вектора.

+ +\param[in] pfft Указатель на структуру `fft_t` алгоритма + быстрого преобразования Фурье.
+ Функция изменит состояние полей структуры `fft_t`, + поэтому структура должна быть очищена перед выходом из + программы для исключения утечек памяти.
+ +\param[in] nfft Размер алгоритма БПФ который будет использован для расчета + секционной свертки с перекрытием.
+ Данный параметр должен быть больше чем минимальное значение + размеров сворачиваемых векторов.
+ Например если `na=10`, а `nb=4`, то параметр `nfft` должен + быть больше 4.
+ Библиотека поддерживает алгоритмы БПФ составной длины + \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). Однако, максимальное + быстродействие достигается при использовании длин равных + степени двойки. + +\param[out] c Указатель на вектор свертки \f$ c = a * b\f$.
+ Размер вектора `[na + nb - 1 x 1]`.
+ Память должна быть выделена.

+ +\return +`RES_OK` если свертка рассчитана успешно.
+ В противном случае \ref ERROR_CODE_GROUP "код ошибки". + + +\note Данная функция наиболее эффективна при вычислении длинных сверток. + +Пример использования функции: + +\include conv_fft_test.c + +Результат работы: +\verbatim + +conv_fft error: 0x00000000 +conv error: 0x00000000 +c[ 0] = -0.00 d[ 0] = 0.00 +c[ 1] = -0.00 d[ 1] = 0.00 +c[ 2] = 1.00 d[ 2] = 1.00 +c[ 3] = 4.00 d[ 3] = 4.00 +c[ 4] = 10.00 d[ 4] = 10.00 +c[ 5] = 20.00 d[ 5] = 20.00 +c[ 6] = 35.00 d[ 6] = 35.00 +c[ 7] = 56.00 d[ 7] = 56.00 +c[ 8] = 77.00 d[ 8] = 77.00 +c[ 9] = 98.00 d[ 9] = 98.00 +c[ 10] = 119.00 d[ 10] = 119.00 +c[ 11] = 140.00 d[ 11] = 140.00 +c[ 12] = 161.00 d[ 12] = 161.00 +c[ 13] = 182.00 d[ 13] = 182.00 +c[ 14] = 190.00 d[ 14] = 190.00 +c[ 15] = 184.00 d[ 15] = 184.00 +c[ 16] = 163.00 d[ 16] = 163.00 +c[ 17] = 126.00 d[ 17] = 126.00 +c[ 18] = 72.00 d[ 18] = 72.00 +\endverbatim + +\author Бахурин Сергей www.dsplib.org +***************************************************************************** */ + + + + + + + + /*! **************************************************************************** \ingroup FILTER_CONV_GROUP \fn int conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb, diff --git a/dspl/src/conv.c b/dspl/src/conv.c index bb1cb87..f47a561 100644 --- a/dspl/src/conv.c +++ b/dspl/src/conv.c @@ -218,18 +218,125 @@ int DSPL_API conv_cmplx(complex_t* a, int na, complex_t* b, free(t); } - return RES_OK; + return RES_OK; } - +/****************************************************************************** +\ingroup FILTER_CONV_GROUP +\fn int conv_fft(double* a, int na, double* b, int nb, + fft_t* pfft, int nfft, double* c) +\brief Real vectors fast linear convolution by using fast Fourier +transform algorithms + +Function convolves two real vectors \f$ c = a * b\f$ length `na` and `nb` +in the frequency domain by using FFT algorithms. This approach provide +high-performance convolution which increases with `na` and `nb` increasing. +The output convolution is a vector `c` with length equal to `na + nb - 1`. + +\param[in] a Pointer to the first vector `a`.
+ Vector size is `[na x 1]`.

+ +\param[in] na Size of the first vector `a`.

+ +\param[in] b Pointer to the second vector `b`.
+ Vector size is `[nb x 1]`.

+ +\param[in] nb Size of the second vector `b`.

+ +\param[in] pfft Pointer to the structure `fft_t`.
+ Function changes `fft_t` structure fields so `fft_t` must + be clear before program returns.

+ +\param[in] nfft FFT size.
+ This parameter set which FFT size will be used + for overlapped frequency domain convolution.
+ FFT size must be more of minimal `na` and `nb` value. + For example if `na = 10`, `nb = 4` then `nfft` parameter must + be more than 4.
+ +\param[out] c Pointer to the convolution output vector \f$ c = a * b\f$.
+ Vector size is `[na + nb - 1 x 1]`.
+ Memory must be allocated.

+ +\return `RES_OK` if convolution is calculated successfully.
+Else \ref ERROR_CODE_GROUP "code error".

+ +Example: +\include conv_fft_test.c + +Program output: + +\verbatim +conv_fft error: 0x00000000 +conv error: 0x00000000 +c[ 0] = -0.00 d[ 0] = 0.00 +c[ 1] = -0.00 d[ 1] = 0.00 +c[ 2] = 1.00 d[ 2] = 1.00 +c[ 3] = 4.00 d[ 3] = 4.00 +c[ 4] = 10.00 d[ 4] = 10.00 +c[ 5] = 20.00 d[ 5] = 20.00 +c[ 6] = 35.00 d[ 6] = 35.00 +c[ 7] = 56.00 d[ 7] = 56.00 +c[ 8] = 77.00 d[ 8] = 77.00 +c[ 9] = 98.00 d[ 9] = 98.00 +c[ 10] = 119.00 d[ 10] = 119.00 +c[ 11] = 140.00 d[ 11] = 140.00 +c[ 12] = 161.00 d[ 12] = 161.00 +c[ 13] = 182.00 d[ 13] = 182.00 +c[ 14] = 190.00 d[ 14] = 190.00 +c[ 15] = 184.00 d[ 15] = 184.00 +c[ 16] = 163.00 d[ 16] = 163.00 +c[ 17] = 126.00 d[ 17] = 126.00 +c[ 18] = 72.00 d[ 18] = 72.00 +\endverbatim + +\author Sergey Bakhurin www.dsplib.org +*******************************************************************************/ +int DSPL_API conv_fft(double* a, int na, double* b, int nb, + fft_t* pfft, int nfft, double* c) +{ + complex_t *pa = NULL, *pb = NULL, *pc = NULL; + int err; + + if(!a || !b || !c || !pfft) + return ERROR_PTR; + if(na<1 || nb < 1) + return ERROR_SIZE; + if(nfft<2) + return ERROR_FFT_SIZE; + + pa = (complex_t*) malloc(na*sizeof(complex_t)); + pb = (complex_t*) malloc(nb*sizeof(complex_t)); + pc = (complex_t*) malloc((na+nb-1)*sizeof(complex_t)); + + re2cmplx(a, na, pa); + re2cmplx(b, nb, pb); + + err = conv_fft_cmplx(pa, na, pb, nb, pfft, nfft, pc); + if(err != RES_OK) + goto exit_label; + + err = cmplx2re(pc, na+nb-1, c, NULL); + +exit_label: + if(pa) free(pa); + if(pb) free(pb); + if(pc) free(pc); + + return err; +} + + + + /****************************************************************************** \ingroup FILTER_CONV_GROUP \fn int conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb, - fft_t* pfft, complex_t* c) + fft_t* pfft, int nfft, complex_t* c) \brief Complex vectors fast linear convolution by using fast Fourier transform algorithms @@ -299,7 +406,7 @@ int DSPL_API conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb, fft_t* pfft, int nfft, complex_t* c) { - int La, Lb, Lc, Nz, n, shift, ind, err; + int La, Lb, Lc, Nz, n, p0, p1, ind, err; complex_t *pa, *pb; complex_t *pt, *pA, *pB, *pC; @@ -312,7 +419,7 @@ int DSPL_API conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb, { La = na; Lb = nb; - pa = a; + pa = a; pb = b; } else @@ -328,77 +435,68 @@ int DSPL_API conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb, if(Nz <= 0) return ERROR_FFT_SIZE; - - + pt = (complex_t*)malloc(nfft*sizeof(complex_t)); pB = (complex_t*)malloc(nfft*sizeof(complex_t)); pA = (complex_t*)malloc(nfft*sizeof(complex_t)); pC = (complex_t*)malloc(nfft*sizeof(complex_t)); - + memset(pt, 0, nfft*sizeof(complex_t)); memcpy(pt+Nz, pb, Lb*sizeof(complex_t)); - + err = fft_cmplx(pt, nfft, pfft, pB); if(err != RES_OK) goto exit_label; - shift = -Lb; + p0 = -Lb; + p1 = p0 + nfft; ind = 0; - while(shift + nfft < La) + while(ind < Lc) { - if(shift < 0) + if(p0 >=0) { - memset(pt, 0, nfft*sizeof(complex_t)); - memcpy(pt-shift, pa, (nfft+shift)*sizeof(complex_t)); - - err = fft_cmplx(pt, nfft, pfft, pA); - if(err != RES_OK) - goto exit_label; + if(p1 < La) + err = fft_cmplx(pa + p0, nfft, pfft, pA); + else + { + memset(pt, 0, nfft*sizeof(complex_t)); + memcpy(pt, pa+p0, (nfft+La-p1)*sizeof(complex_t)); + err = fft_cmplx(pt, nfft, pfft, pA); + } } else { - err = fft_cmplx(pa+shift, nfft, pfft, pA); - if(err != RES_OK) - goto exit_label; + memset(pt, 0, nfft*sizeof(complex_t)); + if(p1 < La) + memcpy(pt - p0, pa, (nfft+p0)*sizeof(complex_t)); + else + memcpy(pt - p0, pa, La * sizeof(complex_t)); + err = fft_cmplx(pt, nfft, pfft, pA); } + + if(err != RES_OK) + goto exit_label; for(n = 0; n < nfft; n++) { RE(pC[n]) = CMRE(pA[n], pB[n]); IM(pC[n]) = CMIM(pA[n], pB[n]); } - - err = ifft_cmplx(pC, nfft, pfft, c+ind); - if(err != RES_OK) - goto exit_label; - - shift += Nz; - ind += Nz; - } - - while(ind <= Lc) - { - memset(pt, 0, nfft*sizeof(complex_t)); - - if(shift>=0) - memcpy(pt, pa + shift, (La-shift)*sizeof(complex_t)); + + + if(ind+nfft < Lc) + err = ifft_cmplx(pC, nfft, pfft, c+ind); else - memcpy(pt+Lb, pa, La*sizeof(complex_t)); - - err = fft_cmplx(pt, nfft, pfft, pA); + { + err = ifft_cmplx(pC, nfft, pfft, pt); + memcpy(c+ind, pt, (Lc-ind)*sizeof(complex_t)); + } if(err != RES_OK) goto exit_label; - for(n = 0; n < nfft; n++) - { - RE(pC[n]) = CMRE(pA[n], pB[n]); - IM(pC[n]) = CMIM(pA[n], pB[n]); - } - - err = ifft_cmplx(pC, nfft, pfft, pt); - memcpy(c+ind, pt, (Lc-ind)*sizeof(complex_t)); - shift += Nz; - ind += Nz; + p0 += Nz; + p1 += Nz; + ind += Nz; } exit_label: @@ -410,167 +508,6 @@ exit_label: return err; } -/*int DSPL_API conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb, - fft_t* pfft, complex_t* c) -{ - complex_t *pa = NULL; - complex_t *pb = NULL; - complex_t *pc = NULL; - complex_t *pA = NULL; - complex_t *pB = NULL; - complex_t *pC = NULL; - - int nfft, nfft2, n, npos, err; - int ma, mb; - complex_t *ta, *tb; - - if(!a || !b || !c) - return ERROR_PTR; - if(na < 1 || nb < 1) - return ERROR_SIZE; - - - if(na > nb) - { - ma = na; - mb = nb; - ta = a; - tb = b; - } - else - { - ma = nb; - mb = na; - ta = b; - tb = a; - } - if(ma > 2*mb) - { - nfft = 4; - n = mb-1; - while(n>>=1) - nfft <<= 1; - nfft2 = nfft >> 1; - - pa = (complex_t*)malloc(nfft * sizeof(complex_t)); - pb = (complex_t*)malloc(nfft * sizeof(complex_t)); - pc = (complex_t*)malloc(nfft * sizeof(complex_t)); - pA = (complex_t*)malloc(nfft * sizeof(complex_t)); - pB = (complex_t*)malloc(nfft * sizeof(complex_t)); - pC = (complex_t*)malloc(nfft * sizeof(complex_t)); - - npos = -nfft2; - memset(pa, 0, nfft*sizeof(complex_t)); - memset(pb, 0, nfft*sizeof(complex_t)); - - memcpy(pa + nfft2, ta, nfft2 * sizeof(complex_t)); - memcpy(pb, tb, mb * sizeof(complex_t)); - - err = fft_cmplx(pa, nfft, pfft, pA); - if(err != RES_OK) - goto exit_label; - - err = fft_cmplx(pb, nfft, pfft, pB); - if(err != RES_OK) - goto exit_label; - - for(n = 0; n < nfft; n++) - { - RE(pC[n]) = CMRE(pA[n], pB[n]); - IM(pC[n]) = CMIM(pA[n], pB[n]); - } - - err = ifft_cmplx(pC, nfft, pfft, pc); - if(err != RES_OK) - goto exit_label; - - memcpy(c, pc+nfft2, nfft2*sizeof(complex_t)); - - npos = 0; - while(npos < ma) - { - if(npos+nfft > ma) - { - memset(pa, 0, nfft * sizeof(complex_t)); - memcpy(pa, ta+npos, (ma - npos) * sizeof(complex_t)); - err = fft_cmplx(pa, nfft, pfft, pA); - - - } - else - err = fft_cmplx(ta+npos, nfft, pfft, pA); - if(err != RES_OK) - goto exit_label; - for(n = 0; n < nfft; n++) - { - RE(pC[n]) = CMRE(pA[n], pB[n]); - IM(pC[n]) = CMIM(pA[n], pB[n]); - } - - err = ifft_cmplx(pC, nfft, pfft, pc); - if(err != RES_OK) - goto exit_label; - if(npos+nfft <= ma+mb-1) - memcpy(c+npos+nfft2, pc+nfft2, - nfft2*sizeof(complex_t)); - else - { - if(ma+mb-1-npos-nfft2 > 0) - { - memcpy(c+npos+nfft2, pc+nfft2,(ma+mb-1-npos-nfft2)*sizeof(complex_t)); - } - } - npos+=nfft2; - } - } - else - { - nfft = 4; - n = ma - 1; - while(n>>=1) - nfft <<= 1; - - pa = (complex_t*)malloc(nfft * sizeof(complex_t)); - pb = (complex_t*)malloc(nfft * sizeof(complex_t)); - pc = (complex_t*)malloc(nfft * sizeof(complex_t)); - pA = (complex_t*)malloc(nfft * sizeof(complex_t)); - pB = (complex_t*)malloc(nfft * sizeof(complex_t)); - pC = (complex_t*)malloc(nfft * sizeof(complex_t)); - - - memset(pa, 0, nfft*sizeof(complex_t)); - memset(pb, 0, nfft*sizeof(complex_t)); - - memcpy(pa, ta, ma * sizeof(complex_t)); - memcpy(pb, tb, mb * sizeof(complex_t)); - - err = fft_cmplx(pa, nfft, pfft, pA); - if(err != RES_OK) - goto exit_label; - err = fft_cmplx(pb, nfft, pfft, pB); - if(err != RES_OK) - goto exit_label; - for(n = 0; n < nfft; n++) - { - RE(pC[n]) = CMRE(pA[n], pB[n]); - IM(pC[n]) = CMIM(pA[n], pB[n]); - } - err = ifft_cmplx(pC, nfft, pfft, pc); - if(err != RES_OK) - goto exit_label; - memcpy(c, pc, (ma+mb-1)*sizeof(complex_t)); - } -exit_label: - if(pa) free(pa); - if(pb) free(pb); - if(pc) free(pc); - if(pA) free(pA); - if(pB) free(pB); - if(pB) free(pC); - - return err; -} -*/ diff --git a/examples/src/conv_fft_test.c b/examples/src/conv_fft_test.c new file mode 100644 index 0000000..0ae3071 --- /dev/null +++ b/examples/src/conv_fft_test.c @@ -0,0 +1,38 @@ +#include +#include +#include +#include "dspl.h" + +#define N 13 +#define M 7 +int main() +{ + void* handle; // DSPL handle + handle = dspl_load(); // Load DSPL function + double a[N], b[M], c[N+M-1], d[N+M-1]; + fft_t pfft; + int n, err; + + linspace(0, N, N, DSPL_PERIODIC, a); + linspace(0, M, M, DSPL_PERIODIC, b); + memset(&pfft, 0, sizeof(fft_t)); + + err = conv_fft(a, N, b, M, &pfft, 16, c); + printf("conv_fft error: 0x%.8x\n", err); + + err = conv(a, N, b, M, d); + printf("conv error: 0x%.8x\n", err); + + // print result + for(n = 0; n < N+M-1; n++) + printf("c[%3d] = %9.2f d[%3d] = %9.2f\n", n, c[n], n, d[n]); + + fft_free(&pfft); // free fft structure memory + dspl_free(handle); // free dspl handle + return 0; +} + + + + + diff --git a/include/dspl.c b/include/dspl.c index fa6fed5..cd1337c 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -52,6 +52,7 @@ p_cmplx2re cmplx2re ; p_concat concat ; p_conv conv ; p_conv_cmplx conv_cmplx ; +p_conv_fft conv_fft ; p_conv_fft_cmplx conv_fft_cmplx ; p_cos_cmplx cos_cmplx ; p_decimate decimate ; @@ -215,6 +216,7 @@ void* dspl_load() LOAD_FUNC(concat); LOAD_FUNC(conv); LOAD_FUNC(conv_cmplx); + LOAD_FUNC(conv_fft); LOAD_FUNC(conv_fft_cmplx); LOAD_FUNC(cos_cmplx); LOAD_FUNC(decimate); diff --git a/include/dspl.h b/include/dspl.h index 190eaf0..c6e4b58 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -344,6 +344,14 @@ DECLARE_FUNC(int, conv_cmplx, complex_t* COMMA int COMMA complex_t*); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, conv_fft, double* a + COMMA int na + COMMA double* b + COMMA int nb + COMMA fft_t* pfft + COMMA int nfft + COMMA double* c); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, conv_fft_cmplx, complex_t* a COMMA int na COMMA complex_t* b diff --git a/release/include/dspl.c b/release/include/dspl.c index fa6fed5..cd1337c 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -52,6 +52,7 @@ p_cmplx2re cmplx2re ; p_concat concat ; p_conv conv ; p_conv_cmplx conv_cmplx ; +p_conv_fft conv_fft ; p_conv_fft_cmplx conv_fft_cmplx ; p_cos_cmplx cos_cmplx ; p_decimate decimate ; @@ -215,6 +216,7 @@ void* dspl_load() LOAD_FUNC(concat); LOAD_FUNC(conv); LOAD_FUNC(conv_cmplx); + LOAD_FUNC(conv_fft); LOAD_FUNC(conv_fft_cmplx); LOAD_FUNC(cos_cmplx); LOAD_FUNC(decimate); diff --git a/release/include/dspl.h b/release/include/dspl.h index 190eaf0..c6e4b58 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -344,6 +344,14 @@ DECLARE_FUNC(int, conv_cmplx, complex_t* COMMA int COMMA complex_t*); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, conv_fft, double* a + COMMA int na + COMMA double* b + COMMA int nb + COMMA fft_t* pfft + COMMA int nfft + COMMA double* c); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, conv_fft_cmplx, complex_t* a COMMA int na COMMA complex_t* b