added kaiser window

pull/6/merge
Dsplib 2019-02-20 23:41:47 +03:00
rodzic 242123b597
commit 1afe24f324
17 zmienionych plików z 517 dodań i 405 usunięć

Wyświetl plik

@ -32,7 +32,6 @@ all:
clean:
$(MAKE) -f Makefile.dspl clean
$(MAKE) -f Makefile.test clean
$(MAKE) -f Makefile.verif clean
$(MAKE) -f Makefile.examples clean

Wyświetl plik

@ -31,7 +31,7 @@
/*! *************************************************************************************************
/*! ****************************************************************************
\ingroup TYPES_GROUP
\def ABSSQR(x)
\brief Макрос возвращает квадрат модуля комплексного числа `x`.
@ -51,12 +51,8 @@
y = ABSSQR(z);
\endcode
Переменная `z = 1-2j`, где `j` - мнимая единица, а переменная `y = 5`.
**************************************************************************************************** */
Переменная `z = 1-2j`, где `j` - мнимая единица, а переменная `y = 5`.
***************************************************************************** */
@ -65,31 +61,30 @@
/*! *************************************************************************************************
/*! ****************************************************************************
\ingroup TYPES_GROUP
\def IM(x)
\brief Макрос определяющий мнимую часть комплексного числа.
Например:
\code
complex_t z;
RE(z) = 1.0;
IM(z) = -2.0;
\endcode
Например:
\code
complex_t z;
RE(z) = 1.0;
IM(z) = -2.0;
\endcode
Переменная `z = 1-2j`, где `j` - мнимая единица.
Переменная `z = 1-2j`, где `j` - мнимая единица.
Аналогично, макрос можно использовать для получения мнимой части комплексного числа:
\code
complex_t z = {3.0, -4.0};
double r;
r = IM(z);
\endcode
В данном примере переменная `z = 3-4i`, а в переменой `r` будет храниться число -4.
**************************************************************************************************** */
Аналогично, макрос можно использовать для получения
мнимой части комплексного числа:
\code
complex_t z = {3.0, -4.0};
double r;
r = IM(z);
\endcode
В данном примере переменная `z = 3-4i`,
а в переменой `r` будет храниться число -4.
***************************************************************************** */
@ -100,31 +95,32 @@
/*! *************************************************************************************************
\ingroup TYPES_GROUP
\def RE(x)
\brief Макрос определяющий реальную часть комплексного числа.
/*! ****************************************************************************
\ingroup TYPES_GROUP
\def RE(x)
\brief Макрос определяющий реальную часть комплексного числа.
Например:
\code
complex_t z;
RE(z) = 1.0;
IM(z) = -2.0;
\endcode
Например:
\code
complex_t z;
RE(z) = 1.0;
IM(z) = -2.0;
\endcode
Переменная `z = 1-2j`, где `j` - мнимая единица.
Переменная `z = 1-2j`, где `j` - мнимая единица.
Аналогично, макрос можно использовать для получения реальной части комплексного числа:
Аналогично, макрос можно использовать для получения
реальной части комплексного числа:
\code
complex_t z = {3.0, -4.0};
double r;
r = RE(z);
\endcode
В данном примере переменная `z = 3-4i`, а в переменой `r` будет храниться число 3.
**************************************************************************************************** */
\code
complex_t z = {3.0, -4.0};
double r;
r = RE(z);
\endcode
В данном примере переменная `z = 3-4i`, а в переменой `r`
будет храниться число 3.
***************************************************************************** */
@ -134,108 +130,104 @@
/*! *************************************************************************************************
\ingroup TYPES_GROUP
\fn int re2cmplx(double* x, int n, complex_t *y)
\brief Преобразование массива вещественных данных в массив комплексных данных.
/*! ****************************************************************************
\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>
Функция заполняет реальные части массива `y` данных соответсвующими значениями
исходного вещественного массива `x`. <BR>
\return
`RES_OK` если преобразование произведено успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки":<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];
Например при выполнении следующего кода
\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
re2cmplx(x, 3, y);
\endcode
\author
Бахурин Сергей.
www.dsplib.org
Значениям `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 Преобразование массива комплексных данных в два массива
вещественных данных, содержащих реальную и мнимую части
исходного массива
/*! ****************************************************************************
\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
Функция заполняет реальные массивы `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

@ -3,159 +3,166 @@
/*! *************************************************************************************************
\ingroup DFT_GROUP
\struct fft_t
\brief Структура данных объекта быстрого преобразования Фурье
Структура хранит указатели на массивы поворотных коэффициентов и массивы промежуточных данных
алгоритма быстрого преобразования Фурье.
Библиотека DSPL использует для БПФ алгоритм для составной длины
/*! ****************************************************************************
\ingroup DFT_GROUP
\struct fft_t
\brief Структура данных объекта быстрого преобразования Фурье
\param n Размер вектора БПФ, для которого выделена память в массивах структуры. <BR>
Парметр `n` должен быть равен целой степени двойки. <BR><BR>
Структура хранит указатели на массивы поворотных коэффициентов
и массивы промежуточных данных алгоритма быстрого преобразования Фурье.
Библиотека DSPL использует для БПФ алгоритм для составной длины
\param w Указатель на вектор поворотных коэффициентов алгоритма БПФ.<BR>
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена и массив поворотных коэффициентов должен
быть заполнен функцией \ref fft_create. <BR><BR>
\param t0 Указатель на вектор промежуточных вычислений алгоритма БПФ.<BR>
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена функцией \ref fft_create. <BR><BR>
\param n Размер вектора БПФ, для которого выделена память
в массивах структуры. <BR>
Парметр `n` должен быть равен целой степени двойки. <BR><BR>
\param t1 Указатель на вектор промежуточных вычислений алгоритма БПФ.<BR>
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена функцией \ref fft_create. <BR><BR>
\param w Указатель на вектор поворотных коэффициентов алгоритма БПФ.<BR>
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена и массив поворотных коэффициентов
должен быть заполнен функцией \ref fft_create. <BR><BR>
Структура заполняется функцией \ref fft_create один раз до использования алгоритма БПФ. <BR>
Указатель на объект данной структуры может быть многократно использован при вызвове функций БПФ.<BR>
Перед выходом из программы выделенную память под поворотные коэффициенты и массивы промежуточных данных
необходимо очистить функцией \ref fft_free. Например:
\param t0 Указатель на вектор промежуточных вычислений алгоритма БПФ.<BR>
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена функцией \ref fft_create. <BR><BR>
\param t1 Указатель на вектор промежуточных вычислений алгоритма БПФ.<BR>
Размер вектора `[n x 1]`. <BR>
Память должна быть выделена функцией \ref fft_create. <BR><BR>
Структура заполняется функцией \ref fft_create один раз
до использования алгоритма БПФ. <BR>
Указатель на объект данной структуры может быть
многократно использован при вызове функций БПФ.<BR>
Перед выходом из программы выделенную память под поворотные
коэффициенты и массивы промежуточных данных
необходимо очистить функцией \ref fft_free. Например:
\code
fft_t pfft; // объявляем объект fft_t
int n = 64; // Размер БПФ
// обнуляем все поля и указатели.
// Данные шаг рекомендуется ввиду того, что некоторые
// при создании переменной не инициализируют ее нулем.
memset(&pfft, 0, sizeof(fft_t));
int err;
//создаем объект для 64-точечного БПФ
err = fft_create(&pfft, n);
// Вызов БПФ функции
// Еще раз вызов БПФ функции
// ....
//очистить память объекта БПФ
fft_free(&pfft);
\endcode
Важно отметить, что если объект `fft_t` был создан для размера БПФ равного `n`,
то он может быть использован только для БПФ размера `n`. <BR>
Также необходимо заметить, что функции БПФ самостоятельно контролируют размер,
и самостоятельно выделяют память объекта БПФ при необходимости.
Так если вызвать любую функцию использующую структуру `fft_t` с заполненными
данными для длины БПФ `k` для расчета БПФ длины `n`,
то массивы стуктуры будут автоматически пересозданы для длины `n`.
\author
Бахурин Сергей.
www.dsplib.org
***************************************************************************** */
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn int fft_create(fft_t *pfft, int n)
\brief Заполнение структуры `fft_t` для алгоритма БПФ
Функция производит выделение памяти и рассчет векторов
поворотных коэффициентов `n`-точечного БПФ для структуры `fft_t`.
\param[in,out] pfft Указатель на структуру `fft_t`. <BR>
Указатель не должен быть `NULL`. <BR><BR>
\param[in] n Размер БПФ \f$n\f$.<BR>
Размер БПФ может быть составным вида
\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.
<BR><BR>
Таким образом алгоритм БПФ поддерживает произвольные
длины, равные целой степени чисел 2,3,5,7,
а также различные их комбинации. <BR><BR>
Так например, при \f$ n = 725760 \f$ структура
будет успешно заполнена, отому что
\f$725760 = 2 \cdot 3 \cdot 4 \cdot 5 \cdot 6 \cdot 7 \cdot 9 \cdot 16 \f$,
т.е. получается как произведение множителей 2,3,5,7.
<BR><BR>
При \f$ n = 172804 = 43201 \cdot 4 \f$ структура также
будет успешно заполнена, потому что простой множитель
входящий в \f$n\f$ не превосходит 46340.
<BR><BR>
Для размера
\f$ n = 13 \cdot 17 \cdot 23 \cdot 13 = 66079 \f$
функция вернет ошибку, поскольку 66079 больше 46340
и не является результатом произведения чисел 2,3,5,7.
<BR><BR>
\return
`RES_OK` если структура заполнена успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".<BR><BR>
\note Некоторые компиляторы при создании структуры не обнуляют ее содержимое.
Поэтому рекомендуется произвести обнуление структуры после ее объявления:
\code
fft_t pfft; // объявляем объект fft_t
int n = 64; // Размер БПФ
\code
fft_t pfft; // объявляем объект fft_t
int n = 64; // Размер БПФ
// обнуляем все поля и указатели.
// Данные шаг рекомендуется ввиду того, что некоторые
// при создании переменной не инициализируют ее нулем.
memset(&pfft, 0, sizeof(fft_t));
int err;
//создаем объект для 64-точечного БПФ
err = fft_create(&pfft, n);
// Вызов БПФ функции
// Еще раз вызов БПФ функции
// ....
//очистить память объекта БПФ
fft_free(&pfft);
\endcode
// обнуляем все поля и указатели.
// Данные шаг рекомендуется ввиду того, что некоторые
// компиляторы при создании переменной не инициализируют ее нулем.
memset(&pfft, 0, sizeof(fft_t));
Важно отметить, что если объект `fft_t` был создан для размера БПФ равного `n`,
то он может быть использован только для БПФ размера `n`. <BR>
int err;
Также необходимо заметить, что функции БПФ самостоятельно контроллируют размер,
и самостоятельно выделяют память объекта БПФ при необходимости. Так если вызвать
любую функцию использующую структуру `fft_t` с заполненными данными для длины
БПФ `k` для расчета БПФ длины `n`, то массивы стуктуры будут автоматически
пересозданы для длины `n`.
\author
Бахурин Сергей.
www.dsplib.org
//создаем объект для 64-точечного БПФ
err = fft_create(&pfft, n);
**************************************************************************************************** */
// .....
//очистить память объекта БПФ
fft_free(&pfft);
\endcode
Перед выходом из программы выделенную в структуре память
необходимо очистить функцией \ref fft_free .<BR><BR>
\note Магия числа 46340 заключается в том, что
\f$\sqrt{2^{31}} = 46340.95\f$.
<BR><BR>
/*! *************************************************************************************************
\ingroup DFT_GROUP
\fn int fft_create(fft_t *pfft, int n)
\brief Заполнение структуры `fft_t` для алгоритма БПФ
Функция производит выделение памяти и рассчет векторов поворотных коэффициентов
`n`-точечного БПФ для структуры `fft_t`.
\param[in,out] pfft Указатель на структуру `fft_t`. <BR>
Указатель не должен быть `NULL`. <BR><BR>
\param[in] n Размер БПФ \f$n\f$.<BR>
Размер БПФ может быть составным вида
\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.<BR><BR>
Таким образом алгоритм БПФ поддерживает произвольный длины, равные целой степени чисел 2,3,5,7,
а таже различные их комбинации. <BR><BR>
Так например, при \f$ n = 725760 \f$ структура будет успешно заполнена,
потому что \f$725760 = 2 \cdot 3 \cdot 4 \cdot 5 \cdot 6 \cdot 7 \cdot 9 \cdot 16 \f$,
т.е. получается как произведение множителей 2,3,5,7.<BR><BR>
При \f$ n = 172804 = 43201 \cdot 4 \f$ структура также будет успешно заполнена,
потому что простой множитель входящий в \f$n\f$ не превосходит 46340.<BR><BR>
Для размера \f$ n = 13 \cdot 17 \cdot 23 \cdot 13 = 66079 \f$ функция вернет ошибку,
поскольку 66079 больше 46340 и не является результатом произведения чисел 2,3,5,7.<BR><BR>
\return
`RES_OK` если структура заполнена успешно. <BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".<BR><BR>
\note Некоторые компиляторы при создании структуры не обнуляют ее содержимое.
Поэтому рекомендуется произвести обнуление структуры после ее объявления:
\code
fft_t pfft; // объявляем объект fft_t
int n = 64; // Размер БПФ
// обнуляем все поля и указатели.
// Данные шаг рекомендуется ввиду того, что некоторые
// компиляторы при создании переменной не инициализируют ее нулем.
memset(&pfft, 0, sizeof(fft_t));
int err;
//создаем объект для 64-точечного БПФ
err = fft_create(&pfft, n);
// .....
//очистить память объекта БПФ
fft_free(&pfft);
\endcode
Перед выходом из программы выделенную в структуре память необходимо очистить функцией \ref fft_free .<BR><BR>
\note Магия числа 46340 заключается в том, что \f$\sqrt{2^{31}} = 46340.95\f$.<BR><BR>
\author
Бахурин Сергей.
www.dsplib.org
**************************************************************************************************** */
\author
Бахурин Сергей.
www.dsplib.org
***************************************************************************** */
@ -165,24 +172,20 @@
/*! *************************************************************************************************
\ingroup DFT_GROUP
\fn void fft_free(fft_t *pfft)
\brief Очистить структуру `fft_t` алгоритма БПФ
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn void fft_free(fft_t *pfft)
\brief Очистить структуру `fft_t` алгоритма БПФ
Функция производит очищение памяти промежуточных данных
и векторов поворотных коэффициентов структуры `fft_t`.
Функция производит очищение памяти промежуточных данных
и векторов поворотных коэффициентов структуры `fft_t`.
\param[in] pfft Указатель на структуру `fft_t`. <BR>
\param[in] pfft Указатель на структуру `fft_t`. <BR>
\author
Бахурин Сергей.
www.dsplib.org
**************************************************************************************************** */
\author
Бахурин Сергей.
www.dsplib.org
***************************************************************************** */
@ -191,38 +194,35 @@
/*! *************************************************************************************************
\ingroup DFT_GROUP
\fn int fft_shift(double* x, int n, double* y)
\brief Перестановка спектральных отсчетов дискретного преобразования Фурье
/*! ****************************************************************************
\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>
Данная функция обрабатывает вещественные входные и выходные вектора и может применяться для перестановки
амплитудного или фазового спектра.
Функция производит
<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>
\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
**************************************************************************************************** */
\author
Бахурин Сергей.
www.dsplib.org
***************************************************************************** */

Wyświetl plik

@ -1,35 +1,84 @@
/*! ****************************************************************************
\ingroup SPEC_MATH_TRANSCEND
\fn int bessel_i0(double* x, int n, double* y)
\brief Модифицированная функция Бесселя первого рода \f$ I_0(x)\f$
Функция рассчитывает значения функции для вещественного вектора `x`,
который должен принимать неотрицательные значения.
<BR>
\param[in] x Указатель на вектор переменной \f$ x \f$.<BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] n Размер входного вектора `x`. <BR><BR>
\param[out] y Указатель на вектор значений функции \f$ I_0(x)\f$.<BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` Расчет произведен успешно.<BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".<BR>
\note Для расчета используется аппроксимация приведенная в статье:<BR>
<a href="https://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision/">
Rational Approximations for the Modified Bessel Function
of the First Kind – I0(x) for Computations with Double Precision<BR>
by PAVEL HOLOBORODKO on NOVEMBER 11, 2015<BR>
</a>
Пример использования функции `bessel_i0`:
\include bessel_i0.c
Данная программа рассчитывает значения функции \f$ I_0(x)\f$ переменной `x`
в интервале \f$[0 \ 3]\f$.
Рассчитанные данные сохраняются в текстовый файл `dat/dat0.txt`
и выводятся на график `img/bessel_i0.png`
\image html bessel_i0.png
Скрипт GNUPLOT для построения графиков из текстовых файлов:
\include bessel_i0.plt
\author
Бахурин Сергей
www.dsplib.org
***************************************************************************** */
/*! ****************************************************************************
\ingroup SPEC_MATH_COMMON_GROUP
\fn int sinc(double* x, int n, double a, double* y)
\brief Функция \f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$
\brief Функция \f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$
Функция рассчитывает занчения значения функции для вещественного вектора `x`.
Функция рассчитывает значения функции для вещественного вектора `x`.
<BR>
\param[in] x Указатель на вектор переменной \f$ x \f$.<BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] x Указатель на вектор переменной \f$ x \f$.<BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] n Размер входного вектора `x`. <BR><BR>
\param[in] n Размер входного вектора `x`. <BR><BR>
\param[in] a Параметр функции
\f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$
\param[in] a Параметр функции
\f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$
\param[out] y Указатель на вектор значений функции.<BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[out] y Указатель на вектор значений функции.<BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` Расчет произведен успешно.<BR>
В противном случае
\ref ERROR_CODE_GROUP "код ошибки".<BR>
`RES_OK` Расчет произведен успешно.<BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".<BR>
\author
Бахурин Сергей
www.dsplib.org
Бахурин Сергей
www.dsplib.org
***************************************************************************** */
@ -78,12 +127,9 @@
Скрипт GNUPLOT для построения графиков из текстовых файлов:
\include sine_int.plt
расчет функции производится на осове
\author
Бахурин Сергей
www.dsplib.org
Бахурин Сергей
www.dsplib.org
***************************************************************************** */

Wyświetl plik

@ -38,6 +38,7 @@
`DSPL_WIN_GAUSSIAN` |Параметрическое окно Гаусса
`DSPL_WIN_HAMMING` |Непараметрическое окно Хемминга
`DSPL_WIN_HANN` |Непараметрическое окно Ханна
`DSPL_WIN_KAISER` |Параметрическое окно Кайзера
`DSPL_WIN_LANCZOS` |Непараметрическое окно Ланкзоса
`DSPL_WIN_NUTTALL` |Непараметрическое окно Натталла
`DSPL_WIN_RECT` |Непараметрическое прямоугольное окно

Wyświetl plik

@ -81,6 +81,7 @@ int win_flat_top (double *w, int n, int win_type);
int win_gaussian (double *w, int n, int win_type, double sigma);
int win_hamming (double *w, int n, int win_type);
int win_hann (double *w, int n, int win_type);
int win_kaiser (double* w, int n, int win_type, double param);
int win_lanczos (double *w, int n, int win_type);
int win_nuttall (double *w, int n, int win_type);
int win_rect (double *w, int n);

Wyświetl plik

@ -33,7 +33,7 @@ Modified Bessel Function of the First Kind – I0(x) [1]
https://www.advanpix.com/2015/11/11/
*******************************************************************************/
int DSPL_API bessel_i0(double *x, int n, double* y)
int DSPL_API bessel_i0(double* x, int n, double* y)
{
double P16[17] = { 1.0000000000000000000000801e+00,
2.4999999999999999999629693e-01,

Wyświetl plik

@ -55,6 +55,8 @@ int window(double* w, int n, int win_type, double param)
return win_hamming(w, n, win_type);
case DSPL_WIN_HANN:
return win_hann(w, n, win_type);
case DSPL_WIN_KAISER:
return win_kaiser(w, n, win_type, param);
case DSPL_WIN_LANCZOS:
return win_lanczos(w, n, win_type);
case DSPL_WIN_NUTTALL:
@ -168,68 +170,6 @@ int win_blackman(double *w, int n, int win_type)
/******************************************************************************
Chebyshev parametric window function
param sets spectrum sidelobes level in dB
ATTENTION! ONLY SYMMETRIC WINDOW
*******************************************************************************/
int win_cheby(double *w, int n, double param)
{
int k, i, m;
double z, dz, sum = 0, wmax=0, r1, x0, chx, chy, in;
if(!w)
return ERROR_PTR;
if(n<2)
return ERROR_SIZE;
if(param <= 0.0)
return ERROR_WIN_PARAM;
r1 = pow(10, param/20);
x0 = cosh((1.0/(double)(n-1)) * acosh(r1));
// check window length even or odd
if(n%2==0)
{
dz = 0.5;
m = n/2-1;
}
else
{
m = (n-1)/2;
dz = 0.0;
}
for(k = 0; k < m+2; k++)
{
z = (double)(k - m) - dz;
sum = 0;
for(i = 1; i <= m; i++)
{
in = (double)i / (double)n;
chx = x0 * cos(M_PI * in);
cheby_poly1(&chx, 1, n-1, &chy);
sum += chy * cos(2.0 * z * M_PI * in);
}
w[k] = r1 + 2.0 * sum;
w[n-1-k] = w[k];
// max value calculation
if(w[k]>wmax)
wmax=w[k];
}
// normalization
for(k=0; k < n; k++)
w[k] /= wmax;
return RES_OK;
}
/******************************************************************************
Blackman - Harris window function
@ -308,6 +248,68 @@ int win_blackman_nuttall(double *w, int n, int win_type)
/******************************************************************************
Chebyshev parametric window function
param sets spectrum sidelobes level in dB
ATTENTION! ONLY SYMMETRIC WINDOW
*******************************************************************************/
int win_cheby(double *w, int n, double param)
{
int k, i, m;
double z, dz, sum = 0, wmax=0, r1, x0, chx, chy, in;
if(!w)
return ERROR_PTR;
if(n<2)
return ERROR_SIZE;
if(param <= 0.0)
return ERROR_WIN_PARAM;
r1 = pow(10, param/20);
x0 = cosh((1.0/(double)(n-1)) * acosh(r1));
// check window length even or odd
if(n%2==0)
{
dz = 0.5;
m = n/2-1;
}
else
{
m = (n-1)/2;
dz = 0.0;
}
for(k = 0; k < m+2; k++)
{
z = (double)(k - m) - dz;
sum = 0;
for(i = 1; i <= m; i++)
{
in = (double)i / (double)n;
chx = x0 * cos(M_PI * in);
cheby_poly1(&chx, 1, n-1, &chy);
sum += chy * cos(2.0 * z * M_PI * in);
}
w[k] = r1 + 2.0 * sum;
w[n-1-k] = w[k];
// max value calculation
if(w[k]>wmax)
wmax=w[k];
}
// normalization
for(k=0; k < n; k++)
w[k] /= wmax;
return RES_OK;
}
/******************************************************************************
@ -486,6 +488,40 @@ int win_hann(double *w, int n, int win_type)
}
/******************************************************************************
Kaiser window function
******************************************************************************/
int win_kaiser(double* w, int n, int win_type, double param)
{
double num, den, x, y;
int i, err;
if(!w)
return ERROR_PTR;
if(n<2)
return ERROR_SIZE;
switch(win_type & DSPL_WIN_SYM_MASK)
{
case DSPL_WIN_SYMMETRIC: x = 1.0/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = 1.0/(double)n; break;
default: return ERROR_WIN_SYM;
}
err = bessel_i0(&param, 1, &den);
if(err != RES_OK)
return err;
for(i = 0; i < n; i++)
{
y = (double)(2*i) / x - 1.0;
y = param * sqrt(1.0 - y*y);
err = bessel_i0(&y, 1, &num);
if(err != RES_OK)
return err;
w[i] = num / den;
}
return err;
}
/******************************************************************************

Wyświetl plik

@ -0,0 +1,11 @@
set grid
set xlabel "x"
set lmargin at screen 0.10
set key left top
set terminal pngcairo size 560,380 enhanced font 'Verdana,8'
set output 'img/bessel_i0.png'
set ylabel "I_0(x)"
set yrange [0:5]
plot 'dat/dat0.txt' with lines

Wyświetl plik

@ -0,0 +1,24 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 50
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double x[N], y[N];
linspace(0.0, 3.0, N, DSPL_SYMMETRIC, x);
bessel_i0(x, N, y);
writetxt(x, y, N, "dat/dat0.txt");
dspl_free(handle); // free dspl handle
return system("gnuplot gnuplot/bessel_i0.plt");
}

Plik binarny nie jest wyświetlany.

Plik binarny nie jest wyświetlany.

Plik binarny nie jest wyświetlany.

Plik binarny nie jest wyświetlany.

Plik binarny nie jest wyświetlany.

Wyświetl plik

@ -173,6 +173,7 @@ typedef struct
#define DSPL_WIN_RECT 0x00010000
#define DSPL_WIN_COS 0x00040000
#define DSPL_WIN_CHEBY 0x00080000
#define DSPL_WIN_KAISER 0x00100000
#define DSPL_FILTER_TYPE_MASK 0x0000000F
#define DSPL_FILTER_LPF 0x00000001

Wyświetl plik

@ -173,6 +173,7 @@ typedef struct
#define DSPL_WIN_RECT 0x00010000
#define DSPL_WIN_COS 0x00040000
#define DSPL_WIN_CHEBY 0x00080000
#define DSPL_WIN_KAISER 0x00100000
#define DSPL_FILTER_TYPE_MASK 0x0000000F
#define DSPL_FILTER_LPF 0x00000001