pull/2/head
Dsplib 2018-05-05 19:51:32 +03:00
rodzic c7dae9b642
commit c79c129fec
11 zmienionych plików z 369 dodań i 26 usunięć

Wyświetl plik

@ -14,31 +14,19 @@
\ref FILTER_ANALYSIS_GROUP <BR>
\ref RESAMPLING_GROUP <BR>
\ref SPEC_MATH_GROUP <BR>
\ref IN_OUT_GROUP <BR>
\ref IN_OUT_GROUP <BR>
<BR>
Библиотека поддерживает работу с вещественными и комплексными типами входных данных. <BR>
Описание испльзуемых типов:<BR>
\ref TYPES_GROUP <BR>
DSPL-2.0 библиотека с открытым исходным кодом, написанная на языке Си.<BR>
DSPL-2.0 использует наиболее эффективные сторонние библиотеки: <BR>
<a href = "http://www.fftw.org/">libfftw3</a> библиотека алгоритмов быстрого преобразования Фурье.<BR>
<a href = "http://www.netlib.org/lapack/">liblapack</a> библиотека алгоритмов линейной алгебры.<BR><BR>
DSPL-2.0 содержит собственные алгоритмы цифровой обработки сигналов, а также реализует интерфейс
для работы со сторонними библиотеками как это показано на структуре: <BR>
\image html dspl_structure.png
В дополнение к этому, совместно с DSPL-2.0, релизована библиотека libplot, которая использует
<a href = "http://ndevilla.free.fr/gnuplot/">ANSI C интерфейс</a> с пакетом построения графиков
<a href = "http://www.gnuplot.info/">GNUPLOT</a>.
Pаспространяется под лицензией <a href = "http://www.gnu.org/licenses/gpl.html">GPL v3</a>
Исходные коды библиотеки доступны на <a href = "https://github.com/Dsplib/dspl-2.0">GitHub</a>.<BR>
Pаспространяется под лицензией <a href = "http://www.gnu.org/licenses/lgpl.html">LGPL v3</a>
Исходные коды библиотеки доступны на <a href = "https://github.com/Dsplib/libdspl-2.0">GitHub</a>.<BR>
Вы также можете внести свой вклад в развитие данной библиотеки. Присоединяйтесь!
*/

Wyświetl plik

@ -1,3 +1,140 @@
/*! *************************************************************************************************
\ingroup TYPES_GROUP
\typedef complex_t
\brief Описание комплексного типа данных.
Комплексный тип данных в библиотеке DSPL-2.0 определен как массив из двух элементов типа `double`.
При этом первый элемент массива определяет реальную часть комплексного числа, а второй - мнимую.
Например:
\code
complex_t z;
z[0] = 1.0;
z[1] = -2.0;
\endcode
Переменная `z = 1-2j`, где `j` - мнимая единица.
Для удобства работы с комплексными числами реализованы специальные макросы:
\ref RE, \ref IM, \ref ABSSQR
**************************************************************************************************** */
/*! *************************************************************************************************
\ingroup TYPES_GROUP
\def ABSSQR(x)
\brief Макрос возвращает квадрат модуля комплексного числа `x`.
Квадрата модуля комплексного числа \f$ x = a + j \cdot b \f$ равен:
\f[
|x|^2 = x \cdot x^* = a^2 + b^2.
\f]
Например:
\code
complex_t z;
double y;
RE(z) = 1.0;
IM(z) = -2.0;
y = ABSSQR(z);
\endcode
Переменная `z = 1-2j`, где `j` - мнимая единица, а переменная `y = 5`.
**************************************************************************************************** */
/*! *************************************************************************************************
\ingroup TYPES_GROUP
\def IM(x)
\brief Макрос определяющий мнимую часть комплексного числа.
Например:
\code
complex_t z;
RE(z) = 1.0;
IM(z) = -2.0;
\endcode
Переменная `z = 1-2j`, где `j` - мнимая единица.
Аналогично, макрос можно использовать для получения мнимой части комплексного числа:
\code
complex_t z = {3.0, -4.0};
double r;
r = IM(z);
\endcode
В данном примере переменная `z = 3-4i`, а в переменой `r` будет храниться число -4.
**************************************************************************************************** */
/*! *************************************************************************************************
\ingroup TYPES_GROUP
\def RE(x)
\brief Макрос определяющий реальную часть комплексного числа.
Например:
\code
complex_t z;
RE(z) = 1.0;
IM(z) = -2.0;
\endcode
Переменная `z = 1-2j`, где `j` - мнимая единица.
Аналогично, макрос можно использовать для получения реальной части комплексного числа:
\code
complex_t z = {3.0, -4.0};
double r;
r = RE(z);
\endcode
В данном примере переменная `z = 3-4i`, а в переменой `r` будет храниться число 3.
**************************************************************************************************** */
/*! *************************************************************************************************
\ingroup TYPES_GROUP

Wyświetl plik

@ -32,12 +32,15 @@
**************************************************************************************************** */
/*! *************************************************************************************************
\ingroup SPEC_MATH_RAND_GEN_GROUP
\fn int randu(double* x, int n);
\brief Генерация вектора равномерно-распределенных в интервале от 0 до 1 псевдослучайных чисел.
Генерация производится при помощи рекурсивного алгоритма L'Ecluyer. Период датчика порядка \f$10^56\f$.<BR>
Генерация производится при помощи рекурсивного алгоритма L'Ecluyer. Период датчика порядка \f$10^{56}\f$.<BR>
\param[in,out] x Указатель на вектор случайных чисел. <BR>

Wyświetl plik

@ -49,8 +49,8 @@ int DSPL_API dft(double* x, int n, complex_t *y)
RE(y[k]) = IM(y[k]) = 0.0;
for(m = 0; m < n; m++)
{
phi = -M_2PI * divn * (double)k * (double)m;
RE(y[k]) += x[m] * cos(phi);
phi = -M_2PI * divn * (double)k * (double)m;
RE(y[k]) += x[m] * cos(phi);
IM(y[k]) += x[m] * sin(phi);
}
}
@ -92,3 +92,44 @@ int DSPL_API dft_cmplx(complex_t* x, int n, complex_t *y)
}
return RES_OK;
}
/**************************************************************************************************
Complex vector DFT
***************************************************************************************************/
int DSPL_API idft_cmplx(complex_t* x, int n, complex_t *y)
{
int k;
int m;
double divn;
double phi;
complex_t e;
if(!x || !y)
return ERROR_PTR;
if(n<1)
return ERROR_SIZE;
divn = 1.0 / (double)n;
for(k = 0; k < n; k++)
{
RE(y[k]) = IM(y[k]) = 0.0;
for(m = 0; m < n; m++)
{
phi = M_2PI * divn * (double)k * (double)m;
RE(e) = cos(phi);
IM(e) = sin(phi);
RE(y[k]) += CMRE(x[m], e);
IM(y[k]) += CMIM(x[m], e);
}
RE(y[k]) /= (double)n;
IM(y[k]) /= (double)n;
}
return RES_OK;
}

Wyświetl plik

@ -29,6 +29,47 @@ int fft_p2(int n);
/**************************************************************************************************
COMPLEX vector IFFT
**************************************************************************************************/
int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
{
int err, k;
double norm;
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));
for(k = 0; k < n; k++)
IM(pfft->t1[k]) = -IM(pfft->t1[k]);
err = fft_dit(pfft, n, y);
if(err!=RES_OK)
return err;
norm = 1.0 / (double)n;
for(k = 0; k < n; k++)
{
RE(y[k]) = RE(y[k])*norm;
IM(y[k]) = -IM(y[k])*norm;
}
return RES_OK;
}
/**************************************************************************************************
Real vector FFT
**************************************************************************************************/
@ -335,6 +376,52 @@ int DSPL_API fft_shift(double* x, int n, double* y)
/**************************************************************************************************
FFT shifting for complex vector
***************************************************************************************************/
int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y)
{
int n2, r;
int k;
complex_t tmp;
complex_t *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++)
{
RE(tmp) = RE(x[k]);
IM(tmp) = IM(x[k]);
RE(y[k]) = RE(x[k+n2]);
IM(y[k]) = IM(x[k+n2]);
RE(y[k+n2]) = RE(tmp);
IM(y[k+n2]) = IM(tmp);
}
}
else
{
n2 = (n-1) >> 1;
buf = (complex_t*) malloc(n2*sizeof(complex_t));
memcpy(buf, x, n2*sizeof(complex_t));
memcpy(y, x+n2, (n2+1)*sizeof(complex_t));
memcpy(y+n2+1, buf, n2*sizeof(complex_t));
free(buf);
}
return RES_OK;
}

Wyświetl plik

@ -92,6 +92,69 @@ exit_label:
/**************************************************************************************************
impulse response of an analog filter H(s)
***************************************************************************************************/
int DSPL_API freqs2time(double* b, double* a, int ord, double fs, int n, fft_t* pfft, double *t, double *h)
{
double *w = NULL;
complex_t *hs = NULL;
complex_t *ht = NULL;
int err, k;
if(!b || !a || !t || !h)
return ERROR_PTR;
if(ord<1)
return ERROR_FILTER_ORD;
if(n<1)
return ERROR_SIZE;
w = (double*)malloc(n*sizeof(double));
hs = (complex_t*)malloc(n*sizeof(complex_t));
err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, w);
if(err != RES_OK)
goto exit_label;
err = freqs(b, a, ord, w, n, hs);
if(err != RES_OK)
goto exit_label;
err = fft_shift_cmplx(hs, n, hs);
if(err != RES_OK)
goto exit_label;
ht = (complex_t*)malloc(n*sizeof(complex_t));
err = ifft_cmplx(hs, n, pfft, ht);
if(err != RES_OK)
{
err = idft_cmplx(hs, n, ht);
if(err != RES_OK)
goto exit_label;
}
for(k = 0; k < n; k++)
{
t[k] = (double)k/fs;
h[k] = RE(ht[k]) * fs;
}
exit_label:
if(w)
free(w);
if(hs)
free(hs);
if(ht)
free(ht);
return err;
}
/**************************************************************************************************
Magnitude, phase response and group delay of an analog filter H(s)

Wyświetl plik

@ -51,6 +51,8 @@ p_cos_cmplx cos_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dmod dmod ;
p_idft_cmplx idft_cmplx ;
p_ifft_cmplx ifft_cmplx ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
p_filter_iir filter_iir ;
@ -60,12 +62,14 @@ p_fft_cmplx fft_cmplx ;
p_fft_create fft_create ;
p_fft_free fft_free ;
p_fft_shift fft_shift ;
p_fft_shift_cmplx fft_shift_cmplx ;
p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_series_dec fourier_series_dec ;
p_fourier_series_rec fourier_series_rec ;
p_freqs freqs ;
p_freqs_resp freqs_resp ;
p_freqs2time freqs2time ;
p_freqz freqz ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
@ -153,6 +157,8 @@ void* dspl_load()
LOAD_FUNC(dft);
LOAD_FUNC(dft_cmplx);
LOAD_FUNC(dmod);
LOAD_FUNC(idft_cmplx);
LOAD_FUNC(ifft_cmplx);
LOAD_FUNC(farrow_lagrange);
LOAD_FUNC(farrow_spline);
LOAD_FUNC(filter_iir);
@ -162,6 +168,7 @@ void* dspl_load()
LOAD_FUNC(fft_create);
LOAD_FUNC(fft_free);
LOAD_FUNC(fft_shift);
LOAD_FUNC(fft_shift_cmplx);
LOAD_FUNC(flipip);
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_series_dec);
@ -169,6 +176,7 @@ void* dspl_load()
LOAD_FUNC(freqz);
LOAD_FUNC(freqs);
LOAD_FUNC(freqs_resp);
LOAD_FUNC(freqs2time);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(linspace);

Wyświetl plik

@ -213,6 +213,8 @@ DECLARE_FUNC(int, cos_cmplx, complex_t* COMMA int COMMA comple
DECLARE_FUNC(int, dft, double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, dft_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(double,dmod, double COMMA double);
DECLARE_FUNC(int, idft_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, ifft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* );
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*);
@ -221,13 +223,15 @@ DECLARE_FUNC(int, fft, double* COMMA int COMMA fft_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, fft_shift, double* COMMA int COMMA double*);
DECLARE_FUNC(int, fft_shift_cmplx, complex_t* COMMA int COMMA complex_t*);
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*);
DECLARE_FUNC(int, fourier_series_rec, double* COMMA complex_t* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqs, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqs_resp, double* COMMA double* COMMA int COMMA double* COMMA int COMMA int COMMA double* COMMA double* COMMA double*);
DECLARE_FUNC(int, freqs2time, double* COMMA double* COMMA int COMMA double COMMA int COMMA fft_t* COMMA double* COMMA double*);
DECLARE_FUNC(int, freqz, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel, double* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel_cmplx, complex_t* COMMA int COMMA int* COMMA int COMMA complex_t*);

Wyświetl plik

@ -51,6 +51,8 @@ p_cos_cmplx cos_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dmod dmod ;
p_idft_cmplx idft_cmplx ;
p_ifft_cmplx ifft_cmplx ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
p_filter_iir filter_iir ;
@ -60,12 +62,14 @@ p_fft_cmplx fft_cmplx ;
p_fft_create fft_create ;
p_fft_free fft_free ;
p_fft_shift fft_shift ;
p_fft_shift_cmplx fft_shift_cmplx ;
p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_series_dec fourier_series_dec ;
p_fourier_series_rec fourier_series_rec ;
p_freqs freqs ;
p_freqs_resp freqs_resp ;
p_freqs2time freqs2time ;
p_freqz freqz ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
@ -153,6 +157,8 @@ void* dspl_load()
LOAD_FUNC(dft);
LOAD_FUNC(dft_cmplx);
LOAD_FUNC(dmod);
LOAD_FUNC(idft_cmplx);
LOAD_FUNC(ifft_cmplx);
LOAD_FUNC(farrow_lagrange);
LOAD_FUNC(farrow_spline);
LOAD_FUNC(filter_iir);
@ -162,6 +168,7 @@ void* dspl_load()
LOAD_FUNC(fft_create);
LOAD_FUNC(fft_free);
LOAD_FUNC(fft_shift);
LOAD_FUNC(fft_shift_cmplx);
LOAD_FUNC(flipip);
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_series_dec);
@ -169,6 +176,7 @@ void* dspl_load()
LOAD_FUNC(freqz);
LOAD_FUNC(freqs);
LOAD_FUNC(freqs_resp);
LOAD_FUNC(freqs2time);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(linspace);

Wyświetl plik

@ -213,6 +213,8 @@ DECLARE_FUNC(int, cos_cmplx, complex_t* COMMA int COMMA comple
DECLARE_FUNC(int, dft, double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, dft_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(double,dmod, double COMMA double);
DECLARE_FUNC(int, idft_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, ifft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* );
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*);
@ -221,13 +223,15 @@ DECLARE_FUNC(int, fft, double* COMMA int COMMA fft_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, fft_shift, double* COMMA int COMMA double*);
DECLARE_FUNC(int, fft_shift_cmplx, complex_t* COMMA int COMMA complex_t*);
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*);
DECLARE_FUNC(int, fourier_series_rec, double* COMMA complex_t* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqs, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, freqs_resp, double* COMMA double* COMMA int COMMA double* COMMA int COMMA int COMMA double* COMMA double* COMMA double*);
DECLARE_FUNC(int, freqs2time, double* COMMA double* COMMA int COMMA double COMMA int COMMA fft_t* COMMA double* COMMA double*);
DECLARE_FUNC(int, freqz, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel, double* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel_cmplx, complex_t* COMMA int COMMA int* COMMA int COMMA complex_t*);

Wyświetl plik

@ -3,7 +3,7 @@
#include <string.h>
#include "dspl.h"
#define ORD 4
#define ORD 3
#define N 1000
@ -13,7 +13,7 @@ int main()
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1];
double Rp = 3.0;
double Rp = 1.0;
double w[N], mag[N], phi[N], tau[N];