added fft and some dox

pull/2/head
Dsplib 2018-03-18 23:11:51 +03:00
rodzic 26c3c73692
commit 8afc35f068
13 zmienionych plików z 613 dodań i 12 usunięć

Wyświetl plik

@ -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)

Wyświetl plik

@ -0,0 +1,109 @@
/*! *************************************************************************************************
\ingroup TYPES_GROUP
\fn int re2cmplx(double* x, int n, complex_t *y)
\brief Преобразование массива вещественных данных в массив комплексных данных.
Функция заполняет реальные части массива `y` данных соответсвующими значениями
исходного вещественного массива `x`. <BR>
\param[in] x Указатель на массв вещественных данных.<BR>
Размер массива `[n x 1]`. <BR><BR>
\param[in] n Размер массивов входных и выходных данных.<BR><BR>
\param[out] y Указатель на адрес массива комплексных данных.<BR>
Размер массива `[n x 1]`. <BR>
Память должна быть выделена. <BR><BR>
\return
`RES_OK` если преобразование произведено успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки":<BR>
Например при выполнении следующего кода
\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`. <BR>
\param[in] x Указатель на массв комплексных данных.<BR>
Размер массива `[n x 1]`. <BR><BR>
\param[in] n Размер массивов входных и выходных данных.<BR><BR>
\param[out] re Указатель на адрес массива реальной части данных.<BR>
Размер массива `[n x 1]`. <BR>
Память должна быть выделена. <BR><BR>
\param[out] im Указатель на адрес массива мнимой части данных.<BR>
Размер массива `[n x 1]`. <BR>
Память должна быть выделена. <BR><BR>
\return
`RES_OK` если преобразование произведено успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки":<BR>
Например при выполнении следующего кода
\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
**************************************************************************************************** */

Wyświetl plik

@ -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 Перестановка спектральных отсчетов дискретного преобразования Фурье
Функция производит
<a href="http://ru.dsplib.org/content/dft_freq/dft_freq.html">
перестановку спектральных отсчетов ДПФ
</a> и переносит нулевую частоту в центр вектора ДПФ. <BR>
Данная функция обрабатывает вещественные входные и выходные вектора и может применяться для перестановки
амплитудного или фазового спектра.
\param[in] x Указатель на исходный вектор ДПФ. <BR>
Размер вектора `[n x 1]`. <BR><BR>
\param[in] n Размер ДПФ \f$n\f$ (размер векторов до и после перестановки).<BR><BR>
\param[out] y Указатель на вектор результата перестановки.<BR>
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` если перестановка произведена успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки":<BR>
\author
Бахурин Сергей.
www.dsplib.org
**************************************************************************************************** */

Wyświetl plik

@ -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

Wyświetl plik

@ -0,0 +1,137 @@
/*!
\ingroup IN_OUT_GROUP
\fn int writebin(void* x, int n, int dtype, char* fn)
\brief Сохранить данные в бинарный файл
Функция сохраняет реальный или комплексный вектор данных
размера `[n x 1]` в бинарный файл `fn`. <BR><BR>
Файл является универсальным для хранения как одномерных,
так и двумерных массивов и имеет следующий формат: <BR><BR>
type 4 байта типа `int`. <BR>
Может принимать значение:<BR>
`DAT_DOUBLE`, если адрес `x` указывает на вектор вещественных чисел; <BR>
`DAT_COMPLEX`, если адрес `x` указывает на вектор комплексных чисел.<BR><BR>
`n` 4 байта типа `int`. <BR>
Количество строк данных. <BR><BR>
`m` 4 байта типа `int`.<BR>
Количество столбцов данных. <BR>
При сохранении вектора всегда равно 1.<BR><BR>
data после идут данные в бинарном виде.<BR>
Размер данных:<BR>
`n * sizeof(double)`, если `dtype==DAT_DOUBLE`;<BR>
`n * sizeof(complex_t)`, если `dtype==DAT_COMPLEX`.<BR><BR>
Файл может быть использован для верификации алгоритмов сторонними пакетами,
такими как GNU Octave или Matlab.<BR><BR>
\param[in] x Указатель на массив данных. <BR>
Размер вектора `[n x 1]`.<BR><BR>
\param[in] n Размер вектора данных.<BR><BR>
\param[in] dtype Тип данных.<BR>
Может принимать значения: `DAT_DOUBLE` или `DAT_COMPLEX`.<BR><BR>
\param[in] fn Имя файла.<BR><BR>
\return
`RES_OK` Файл сохранен успешно.<BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки":<BR>
\note
Данная функция производит запись в файл без потери точности,
поэтому рекомендуется использовать ее для верификации данных DSPL.<BR><BR>
Функция для чтения бинарного файла в 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`. <BR>
Файл имеет следующий формат<BR>
\verbatim
x[0] y[0]
x[1] y[1]
...
x[n-1] y[n-1]
\endverbatim
Файл может быть использован для построения графика сторонней программой.<BR><BR>
\param[in] x Указатель на первый вектор. <BR>
Размер вектора `[n x 1]`.<BR><BR>
\param[in] y Указатель на второй вектор. <BR>
Размер вектора `[n x 1]`. <BR>
Может быть `NULL`. <BR>
Файл будет содержать только один столбец соответствующий
вектору `x` если `y == NULL`.<BR><BR>
\param[in] n Размер входных векторов.<BR><BR>
\param[in] fn Имя файла.<BR><BR>
\return
`RES_OK` Файл сохранен успешно.<BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки":<BR>
\note
Данная функция производит округление данных при записи в файл.
Поэтому не рекомендуется использовать ее для верификации данных DSPL.
\author
Бахурин Сергей.
www.dsplib.org
*/

Wyświetl plik

@ -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<<m;
for( r = 0; r < s; r++)
{
phi = -M_2PI *(double)r / (double)(2*s);
RE(pfft->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;
}

Wyświetl plik

@ -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);

Wyświetl plik

@ -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*);

Wyświetl plik

@ -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);

Wyświetl plik

@ -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*);

Plik binarny nie jest wyświetlany.

Wyświetl plik

@ -14,6 +14,7 @@ int main()
//complex_t z[N];
for(int k = 0; k < N; k++)
{
RE(x[k]) = (double)k;

Wyświetl plik

@ -0,0 +1,36 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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;
}