diff --git a/Makefile.test b/Makefile.test index 334fecb..2deffda 100644 --- a/Makefile.test +++ b/Makefile.test @@ -1,8 +1,8 @@ -PRJ_DIR = test - +PRJ_DIR = test SRC_DIR = $(PRJ_DIR)/src BIN_DIR = $(PRJ_DIR)/bin + include Makefile.dirs DSPL_C_FILE = $(INC_DIR)/dspl.c @@ -11,12 +11,11 @@ DSPL_O_FILE = $(PRJ_DIR)/obj/dspl.o SRC_FILES = $(wildcard $(SRC_DIR)/*.c) BIN_FILES = $(addprefix $(BIN_DIR)/,$(notdir $(SRC_FILES:.c=.exe))) - CFLAGS = -Wall -O3 -I$(INC_DIR) -D$(DEF_OS) all: $(BIN_FILES) -$(BIN_DIR)/%.exe: $(SRC_DIR)/%.c $(DSPL_O_FILE) +$(BIN_DIR)/%.exe:$(SRC_DIR)/%.c $(DSPL_O_FILE) $(CC) $(CFLAGS) $< $(DSPL_O_FILE) -o $@ $(LFLAGS) $(DSPL_O_FILE):$(DSPL_C_FILE) diff --git a/dspl/dox/ru/complex.dox b/dspl/dox/ru/complex.dox new file mode 100644 index 0000000..af1ec83 --- /dev/null +++ b/dspl/dox/ru/complex.dox @@ -0,0 +1,109 @@ + +/*! ************************************************************************************************* + \ingroup TYPES_GROUP + \fn int re2cmplx(double* x, int n, complex_t *y) + \brief Преобразование массива вещественных данных в массив комплексных данных. + + Функция заполняет реальные части массива `y` данных соответсвующими значениями + исходного вещественного массива `x`.
+ + + \param[in] x Указатель на массв вещественных данных.
+ Размер массива `[n x 1]`.

+ + \param[in] n Размер массивов входных и выходных данных.

+ + \param[out] y Указатель на адрес массива комплексных данных.
+ Размер массива `[n x 1]`.
+ Память должна быть выделена.

+ + + \return + `RES_OK` если преобразование произведено успешно.
+ В противном случае \ref ERROR_CODE_GROUP "код ошибки":
+ + + + Например при выполнении следующего кода + \code{.cpp} + double x[3] = {1.0, 2.0, 3.0}; + complex_t y[3]; + + re2cmplx(x, 3, y); + \endcode + + Значениям `y` будут присвоены значения: + + \verbatim + y[0] = 1+0j; + y[1] = 2+0j; + y[2] = 3+0j. + \endverbatim + + \author + Бахурин Сергей. + www.dsplib.org + + +**************************************************************************************************** */ + + + + + + +/*! ************************************************************************************************* + \ingroup TYPES_GROUP + \fn int cmplx2re(complex_t* x, int n, double *re, double *im) + \brief Преобразование массива комплексных данных в два массива + вещественных данных, содержащих реальную и мнимую части + исходного массива + + Функция заполняет реальные массивы `re` и `im` соответсвующими значениями + ральной и мнимой частей исходного комплексного массива `x`.
+ + + \param[in] x Указатель на массв комплексных данных.
+ Размер массива `[n x 1]`.

+ + \param[in] n Размер массивов входных и выходных данных.

+ + \param[out] re Указатель на адрес массива реальной части данных.
+ Размер массива `[n x 1]`.
+ Память должна быть выделена.

+ + \param[out] im Указатель на адрес массива мнимой части данных.
+ Размер массива `[n x 1]`.
+ Память должна быть выделена.

+ + \return + `RES_OK` если преобразование произведено успешно.
+ В противном случае \ref ERROR_CODE_GROUP "код ошибки":
+ + + + Например при выполнении следующего кода + \code{.cpp} + complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}}; + double re[3], im[3]; + + cmplx2re(x, 3, re, im); + \endcode + + Элементам массивов `re` и `im` будут присвоены значения: + + + \verbatim + re[0] = 1.0; im[0] = 2.0; + re[1] = 3.0; im[1] = 4.0; + re[2] = 5.0; im[2] = 6.0; + \endverbatim + + \author + Бахурин Сергей. + www.dsplib.org + + +**************************************************************************************************** */ + + diff --git a/dspl/dox/ru/fft.dox b/dspl/dox/ru/fft.dox index 1e5a0fe..e099c3f 100644 --- a/dspl/dox/ru/fft.dox +++ b/dspl/dox/ru/fft.dox @@ -1,5 +1,37 @@ -/************************************************************************************************** -FFT shifting -int fft_shift(double* x, int n, double* y) -***************************************************************************************************/ + +/*! ************************************************************************************************* + \ingroup DFT_GROUP + \fn int fft_shift(double* x, int n, double* y) + \brief Перестановка спектральных отсчетов дискретного преобразования Фурье + + Функция производит + + перестановку спектральных отсчетов ДПФ + и переносит нулевую частоту в центр вектора ДПФ.
+ Данная функция обрабатывает вещественные входные и выходные вектора и может применяться для перестановки + амплитудного или фазового спектра. + + + + + \param[in] x Указатель на исходный вектор ДПФ.
+ Размер вектора `[n x 1]`.

+ + \param[in] n Размер ДПФ \f$n\f$ (размер векторов до и после перестановки).

+ + \param[out] y Указатель на вектор результата перестановки.
+ Размер вектора `[n x 1]`.
+ Память должна быть выделена.

+ + + \return + `RES_OK` если перестановка произведена успешно.
+ В противном случае \ref ERROR_CODE_GROUP "код ошибки":
+ + \author + Бахурин Сергей. + www.dsplib.org + + +**************************************************************************************************** */ diff --git a/dspl/dox/ru/fillarray.dox b/dspl/dox/ru/fillarray.dox index 48e2888..9f80f3f 100644 --- a/dspl/dox/ru/fillarray.dox +++ b/dspl/dox/ru/fillarray.dox @@ -38,7 +38,7 @@ Пример 1. Периодическое заполнение. \code double x[5]; - dspl_linspace(0, 5, 5, DSPL_PERIODIC, x); + linspace(0, 5, 5, DSPL_PERIODIC, x); \endcode В массиве `x` будут лежать значения: \code @@ -48,7 +48,7 @@ Пример 2. Симметричное заполнение. \code double x[5]; - dspl_linspace(0, 5, 5, DSPL_SYMMETRIC, x); + linspace(0, 5, 5, DSPL_SYMMETRIC, x); \endcode В массиве `x` будут лежать значения: \code @@ -108,7 +108,7 @@ Пример 1. Периодическое заполнение. \code double x[5]; - dspl_logspace(-2, 3, 5, DSPL_PERIODIC, x); + logspace(-2, 3, 5, DSPL_PERIODIC, x); \endcode В массиве `x` будут лежать значения: \code @@ -118,7 +118,7 @@ Пример 2. Симметричное заполнение. \code double x[5]; - dspl_logspace(-2, 3, 5, DSPL_SYMMETRIC, x); + logspace(-2, 3, 5, DSPL_SYMMETRIC, x); \endcode В массиве `x` будут лежать значения: \code diff --git a/dspl/dox/ru/inout.dox b/dspl/dox/ru/inout.dox new file mode 100644 index 0000000..33937b0 --- /dev/null +++ b/dspl/dox/ru/inout.dox @@ -0,0 +1,137 @@ +/*! + \ingroup IN_OUT_GROUP + \fn int writebin(void* x, int n, int dtype, char* fn) + \brief Сохранить данные в бинарный файл + + Функция сохраняет реальный или комплексный вектор данных + размера `[n x 1]` в бинарный файл `fn`.

+ + Файл является универсальным для хранения как одномерных, + так и двумерных массивов и имеет следующий формат:

+ + type 4 байта типа `int`.
+ Может принимать значение:
+ `DAT_DOUBLE`, если адрес `x` указывает на вектор вещественных чисел;
+ `DAT_COMPLEX`, если адрес `x` указывает на вектор комплексных чисел.

+ + `n` 4 байта типа `int`.
+ Количество строк данных.

+ + `m` 4 байта типа `int`.
+ Количество столбцов данных.
+ При сохранении вектора всегда равно 1.

+ + data после идут данные в бинарном виде.
+ Размер данных:
+ `n * sizeof(double)`, если `dtype==DAT_DOUBLE`;
+ `n * sizeof(complex_t)`, если `dtype==DAT_COMPLEX`.

+ + Файл может быть использован для верификации алгоритмов сторонними пакетами, + такими как GNU Octave или Matlab.

+ + \param[in] x Указатель на массив данных.
+ Размер вектора `[n x 1]`.

+ + \param[in] n Размер вектора данных.

+ + \param[in] dtype Тип данных.
+ Может принимать значения: `DAT_DOUBLE` или `DAT_COMPLEX`.

+ + \param[in] fn Имя файла.

+ + \return + `RES_OK` Файл сохранен успешно.
+ В противном случае \ref ERROR_CODE_GROUP "код ошибки":
+ + \note + Данная функция производит запись в файл без потери точности, + поэтому рекомендуется использовать ее для верификации данных DSPL.

+ + Функция для чтения бинарного файла в GNU Octave и Matlab: + + \code + function [dat, n, m] = readbin(fn) + + fid = fopen(fn); + if(~fid) + error('cannot to open file'); + end + type = fread(fid, 1, 'int32'); + n = fread(fid, 1, 'int32'); + m = fread(fid, 1, 'int32'); + + if(type==0) + dat = fread(fid, [n*m, 1], 'double'); + end + + if(type==1) + y = fread(fid, [n*m*2, 1], 'double'); + dat = y(1:2:end) + 1i * y(2:2:end); + end + + dat = reshape(dat, n, m); + + fclose(fid); + end + \endcode + + + \author + Бахурин Сергей. + www.dsplib.org + + +*/ + + + + + + + +/*! + \ingroup IN_OUT_GROUP + \fn int writetxt(double* x, double *y, int n, char* fn) + \brief Сохранить вещественные данные в текстовый файл + + Функция сохраняет вещественные данные в текстовый файл `fn`.
+ + Файл имеет следующий формат
+ + \verbatim + x[0] y[0] + x[1] y[1] + ... + x[n-1] y[n-1] + \endverbatim + + Файл может быть использован для построения графика сторонней программой.

+ + \param[in] x Указатель на первый вектор.
+ Размер вектора `[n x 1]`.

+ + \param[in] y Указатель на второй вектор.
+ Размер вектора `[n x 1]`.
+ Может быть `NULL`.
+ Файл будет содержать только один столбец соответствующий + вектору `x` если `y == NULL`.

+ + \param[in] n Размер входных векторов.

+ + \param[in] fn Имя файла.

+ + + \return + `RES_OK` Файл сохранен успешно.
+ В противном случае \ref ERROR_CODE_GROUP "код ошибки":
+ + \note + Данная функция производит округление данных при записи в файл. + Поэтому не рекомендуется использовать ее для верификации данных DSPL. + + \author + Бахурин Сергей. + www.dsplib.org + + +*/ diff --git a/dspl/src/fft.c b/dspl/src/fft.c index 820b25c..3888ae7 100644 --- a/dspl/src/fft.c +++ b/dspl/src/fft.c @@ -23,10 +23,60 @@ #include "dspl.h" int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2); +int fft_dit(fft_t *pfft, int n, complex_t* y); +void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n, complex_t *y0, complex_t *y1); int fft_p2(int n); +/************************************************************************************************** +Real vector FFT +**************************************************************************************************/ +int DSPL_API fft(double *x, int n, fft_t* pfft, complex_t* y) +{ + int err; + + if(!x || !pfft || !y) + return ERROR_PTR; + if(n<1) + return ERROR_SIZE; + + + err = fft_create(pfft, n); + if(err != RES_OK) + return err; + + re2cmplx(x, n, pfft->t1); + + return fft_dit(pfft, n, y); +} + + + + +/************************************************************************************************** +COMPLEX vector FFT +**************************************************************************************************/ +int DSPL_API fft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y) +{ + int err; + + if(!x || !pfft || !y) + return ERROR_PTR; + if(n<1) + return ERROR_SIZE; + + + err = fft_create(pfft, n); + if(err != RES_OK) + return err; + + memcpy(pfft->t1, x, n*sizeof(complex_t)); + + return fft_dit(pfft, n, y); +} + + /************************************************************************************************** FFT bit reverse @@ -77,6 +127,150 @@ int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2) +/************************************************************************************************** +FFT create +**************************************************************************************************/ +int DSPL_API fft_create(fft_t *pfft, int n) +{ + int p2, k,r,m,addr,s; + double phi; + + p2 = fft_p2(n); + if(p2 < 1) + return ERROR_FFT_SIZE; + if(n < pfft->n+1) + return RES_OK; + + pfft->n = n; + pfft->p2 = p2; + + pfft->w = pfft->w ? (complex_t*) realloc(pfft->w, n*sizeof(complex_t)): + (complex_t*) malloc( n*sizeof(complex_t)); + + + pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)): + (complex_t*) malloc( n*sizeof(complex_t)); + + + pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)): + (complex_t*) malloc( n*sizeof(complex_t)); + + m = 0; + addr = 0; + + for(k = 0; k < p2; k++) + { + s = 1<w[addr+r]) = cos(phi); + IM(pfft->w[addr+r]) = sin(phi); + } + addr+=s; + m++; + } + return RES_OK; +} + + + +/************************************************************************************************** +FFT decimation in time +**************************************************************************************************/ +int fft_dit(fft_t *pfft, int n, complex_t* y) +{ + int k,s,m,waddr, dm,p2; + complex_t *t, *t0, *t1, *w; + int err; + + p2 = fft_p2(n); + if(p2<0) + return ERROR_FFT_SIZE; + + t0 = pfft->t0; + t1 = pfft->t1; + w = pfft->w; + + s = n>>1; + m = 1; + waddr = 0; + + err = fft_bit_reverse(t1, t0, n, p2); + if(err!= RES_OK) + return err; + + while(s) + { + dm = m<<1; + if(s > 1) + { + for(k = 0; k < n; k+=dm) + { + fft_dit_krn(t0+k, t0+k+m, w+waddr, m, t1+k, t1+k+m); + t = t1; + t1 = t0; + t0 = t; + waddr+=m; + m <<= 1; + } + } + else + { + fft_dit_krn(t0, t0+m, w+waddr, m, y, y+m); + } + s >>= 1; + } + + + return RES_OK; +} + + + + + +/************************************************************************************************** +FFT decimation in time kernel +**************************************************************************************************/ +void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n, complex_t *y0, complex_t *y1) +{ + int k; + complex_t mul; + for(k = 0; k < n; k++) + { + RE(mul) = CMRE(x1[k], w[k]); + IM(mul) = CMIM(x1[k], w[k]); + + RE(y0[k]) = RE(x0[k]) + RE(mul); + IM(y0[k]) = IM(x0[k]) + IM(mul); + + RE(y1[k]) = RE(x0[k]) - RE(mul); + IM(y1[k]) = IM(x0[k]) - IM(mul); + } + +} + + + + + +/************************************************************************************************** +FFT free +**************************************************************************************************/ +void DSPL_API fft_free(fft_t *pfft) +{ + if(!pfft) + return; + if(pfft->w) + free(pfft->w); + if(pfft->t0) + free(pfft->t0); + if(pfft->t1) + free(pfft->t1); +} + + @@ -93,3 +287,56 @@ int fft_p2(int n) return -1; return p2; } + + + + + + + +/************************************************************************************************** +FFT shifting +***************************************************************************************************/ +int DSPL_API fft_shift(double* x, int n, double* y) +{ + int n2, r; + int k; + double tmp; + double *buf; + + if(!x || !y) + return ERROR_PTR; + + if(n<1) + return ERROR_SIZE; + + r = n%2; + if(!r) + { + n2 = n>>1; + for(k = 0; k < n2; k++) + { + tmp = x[k]; + y[k] = x[k+n2]; + y[k+n2] = tmp; + } + } + else + { + n2 = (n-1) >> 1; + buf = (double*) malloc(n2*sizeof(double)); + memcpy(buf, x, n2*sizeof(double)); + memcpy(y, x+n2, (n2+1)*sizeof(double)); + memcpy(y+n2+1, buf, n2*sizeof(double)); + free(buf); + } + return RES_OK; +} + + + + + + + + diff --git a/include/dspl.c b/include/dspl.c index 38e1d06..ec12f5b 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -35,6 +35,8 @@ #ifndef BUILD_LIB + + p_acos_cmplx acos_cmplx ; p_asin_cmplx asin_cmplx ; p_cheby_poly1 cheby_poly1 ; @@ -50,6 +52,11 @@ p_dmod dmod ; p_farrow_lagrange farrow_lagrange ; p_farrow_spline farrow_spline ; p_filter_iir filter_iir ; +p_fft fft ; +p_fft_cmplx fft_cmplx ; +p_fft_create fft_create ; +p_fft_free fft_free ; +p_fft_shift fft_shift ; p_flipip flipip ; p_flipip_cmplx flipip_cmplx ; p_fourier_series_dec fourier_series_dec ; @@ -137,6 +144,11 @@ void* dspl_load() LOAD_FUNC(farrow_lagrange); LOAD_FUNC(farrow_spline); LOAD_FUNC(filter_iir); + LOAD_FUNC(fft); + LOAD_FUNC(fft_cmplx); + LOAD_FUNC(fft_create); + LOAD_FUNC(fft_free); + LOAD_FUNC(fft_shift); LOAD_FUNC(flipip); LOAD_FUNC(flipip_cmplx); LOAD_FUNC(fourier_series_dec); diff --git a/include/dspl.h b/include/dspl.h index adf2c99..8805244 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -82,6 +82,7 @@ typedef struct #define ERROR_ELLIP_MODULE 0x05121315 /* F 0x06xxxxxx*/ #define ERROR_FFT_CREATE 0x06060318 +#define ERROR_FFT_SIZE 0x06062021 #define ERROR_FILTER_A0 0x06100100 #define ERROR_FILTER_ORD 0x06101518 #define ERROR_FILTER_RP 0x06101816 @@ -206,6 +207,13 @@ DECLARE_FUNC(double,dmod, double COMMA double); DECLARE_FUNC(int, farrow_lagrange, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*); DECLARE_FUNC(int, farrow_spline, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*); DECLARE_FUNC(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*); + + +DECLARE_FUNC(int, fft, double* COMMA int COMMA fft_t* COMMA complex_t* ); +DECLARE_FUNC(int, fft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* ); +DECLARE_FUNC(int, fft_create, fft_t* COMMA int); +DECLARE_FUNC(void, fft_free, fft_t*); +DECLARE_FUNC(int, fft_shift, double* COMMA int n COMMA double*); DECLARE_FUNC(int, flipip, double* COMMA int); DECLARE_FUNC(int, flipip_cmplx, complex_t* COMMA int); DECLARE_FUNC(int, fourier_series_dec, double* COMMA double* COMMA int COMMA double COMMA int COMMA double* COMMA complex_t*); diff --git a/release/include/dspl.c b/release/include/dspl.c index 38e1d06..ec12f5b 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -35,6 +35,8 @@ #ifndef BUILD_LIB + + p_acos_cmplx acos_cmplx ; p_asin_cmplx asin_cmplx ; p_cheby_poly1 cheby_poly1 ; @@ -50,6 +52,11 @@ p_dmod dmod ; p_farrow_lagrange farrow_lagrange ; p_farrow_spline farrow_spline ; p_filter_iir filter_iir ; +p_fft fft ; +p_fft_cmplx fft_cmplx ; +p_fft_create fft_create ; +p_fft_free fft_free ; +p_fft_shift fft_shift ; p_flipip flipip ; p_flipip_cmplx flipip_cmplx ; p_fourier_series_dec fourier_series_dec ; @@ -137,6 +144,11 @@ void* dspl_load() LOAD_FUNC(farrow_lagrange); LOAD_FUNC(farrow_spline); LOAD_FUNC(filter_iir); + LOAD_FUNC(fft); + LOAD_FUNC(fft_cmplx); + LOAD_FUNC(fft_create); + LOAD_FUNC(fft_free); + LOAD_FUNC(fft_shift); LOAD_FUNC(flipip); LOAD_FUNC(flipip_cmplx); LOAD_FUNC(fourier_series_dec); diff --git a/release/include/dspl.h b/release/include/dspl.h index adf2c99..8805244 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -82,6 +82,7 @@ typedef struct #define ERROR_ELLIP_MODULE 0x05121315 /* F 0x06xxxxxx*/ #define ERROR_FFT_CREATE 0x06060318 +#define ERROR_FFT_SIZE 0x06062021 #define ERROR_FILTER_A0 0x06100100 #define ERROR_FILTER_ORD 0x06101518 #define ERROR_FILTER_RP 0x06101816 @@ -206,6 +207,13 @@ DECLARE_FUNC(double,dmod, double COMMA double); DECLARE_FUNC(int, farrow_lagrange, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*); DECLARE_FUNC(int, farrow_spline, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*); DECLARE_FUNC(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*); + + +DECLARE_FUNC(int, fft, double* COMMA int COMMA fft_t* COMMA complex_t* ); +DECLARE_FUNC(int, fft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* ); +DECLARE_FUNC(int, fft_create, fft_t* COMMA int); +DECLARE_FUNC(void, fft_free, fft_t*); +DECLARE_FUNC(int, fft_shift, double* COMMA int n COMMA double*); DECLARE_FUNC(int, flipip, double* COMMA int); DECLARE_FUNC(int, flipip_cmplx, complex_t* COMMA int); DECLARE_FUNC(int, fourier_series_dec, double* COMMA double* COMMA int COMMA double COMMA int COMMA double* COMMA complex_t*); diff --git a/test/bin/dft_test.exe b/test/bin/dft_test.exe index de5dff4..3f35fe3 100755 Binary files a/test/bin/dft_test.exe and b/test/bin/dft_test.exe differ diff --git a/test/src/dft_test.c b/test/src/dft_test.c index f6feb10..6825f93 100644 --- a/test/src/dft_test.c +++ b/test/src/dft_test.c @@ -14,6 +14,7 @@ int main() //complex_t z[N]; + for(int k = 0; k < N; k++) { RE(x[k]) = (double)k; diff --git a/test/src/fft_test.c b/test/src/fft_test.c new file mode 100644 index 0000000..597cd4c --- /dev/null +++ b/test/src/fft_test.c @@ -0,0 +1,36 @@ +#include +#include +#include +#include "dspl.h" + +#define N 16 +int main() +{ + void* handle; + handle = dspl_load(); + + complex_t x[N]; + complex_t y[N]; + fft_t pfft; + + for(int k = 0; k < N; k++) + { + RE(x[k]) = (double)k; + IM(x[k]) = 0.0; + } + memset(&pfft, 0, sizeof(fft_t)); + + fft_create(&pfft, N); + + fft_cmplx(x, N, &pfft, y); + + //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); + + dspl_free(handle); + return 0; +} + +