diff --git a/dspl/dox/ru/conv.dox b/dspl/dox/ru/conv.dox index 881a12f..78283d7 100644 --- a/dspl/dox/ru/conv.dox +++ b/dspl/dox/ru/conv.dox @@ -163,6 +163,20 @@ cc[5] = 0.0+24.0j поэтому структура должна быть очищена перед выходом из программы для исключения утечек памяти.
+\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]`.
Память должна быть выделена.

diff --git a/dspl/src/conv.c b/dspl/src/conv.c index 8467be5..bb1cb87 100644 --- a/dspl/src/conv.c +++ b/dspl/src/conv.c @@ -251,6 +251,13 @@ The output convolution is a vector `c` with length equal to `na + nb - 1`. \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]`.
@@ -288,7 +295,122 @@ c[ 18] = -37.00 +484.00j d[ 18] = -37.00 +484.00j \author Sergey Bakhurin www.dsplib.org *******************************************************************************/ -int DSPL_API conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb, +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; + complex_t *pa, *pb; + complex_t *pt, *pA, *pB, *pC; + + if(!a || !b || !c) + return ERROR_PTR; + if(na < 1 || nb < 1) + return ERROR_SIZE; + + if(na >= nb) + { + La = na; + Lb = nb; + pa = a; + pb = b; + } + else + { + La = nb; + pa = b; + Lb = na; + pb = a; + } + + Lc = La + Lb - 1; + Nz = nfft - Lb; + + 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; + ind = 0; + while(shift + nfft < La) + { + if(shift < 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; + } + else + { + err = fft_cmplx(pa+shift, 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)); + else + memcpy(pt+Lb, 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, pt); + memcpy(c+ind, pt, (Lc-ind)*sizeof(complex_t)); + shift += Nz; + ind += Nz; + } + +exit_label: + if(pt) free(pt); + if(pB) free(pB); + if(pA) free(pA); + if(pC) free(pC); + + 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; @@ -448,7 +570,7 @@ exit_label: return err; } - +*/ diff --git a/examples/src/conv_fft_cmplx_test.c b/examples/src/conv_fft_cmplx_test.c index 16601b9..3e9f810 100644 --- a/examples/src/conv_fft_cmplx_test.c +++ b/examples/src/conv_fft_cmplx_test.c @@ -17,7 +17,7 @@ int main() linspace(0, 2*M, 2*M, DSPL_PERIODIC, (double*)b); memset(&pfft, 0, sizeof(fft_t)); - conv_fft_cmplx(a, N, b, M, &pfft, c); + conv_fft_cmplx(a, N, b, M, &pfft, 8, c); conv_cmplx(a, N, b, M, d); // print result diff --git a/include/dspl.h b/include/dspl.h index 07ae6a3..4e4f814 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -349,6 +349,7 @@ DECLARE_FUNC(int, conv_fft_cmplx, complex_t* a COMMA complex_t* b COMMA int nb COMMA fft_t* pfft + COMMA int nfft COMMA complex_t* c); //------------------------------------------------------------------------------ DECLARE_FUNC(int, cos_cmplx, complex_t* diff --git a/release/include/dspl.h b/release/include/dspl.h index 07ae6a3..4e4f814 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -349,6 +349,7 @@ DECLARE_FUNC(int, conv_fft_cmplx, complex_t* a COMMA complex_t* b COMMA int nb COMMA fft_t* pfft + COMMA int nfft COMMA complex_t* c); //------------------------------------------------------------------------------ DECLARE_FUNC(int, cos_cmplx, complex_t*