kopia lustrzana https://github.com/Dsplib/libdspl-2.0
added goertzel alg, fft bit reverse and power of 2. Also added polyval functions and some docs
Dsplib
rodzic
ff3019dd5a
commit
3daecefae2
|
|
@ -781,7 +781,12 @@ WARN_LOGFILE =
|
|||
# spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
|
||||
# Note: If this tag is empty the current directory is searched.
|
||||
|
||||
INPUT = ../
|
||||
INPUT = ru \
|
||||
../dspl/dox/ru \
|
||||
../dspl/src \
|
||||
../include \
|
||||
../test/dox/ru \
|
||||
../test/src
|
||||
|
||||
# This tag can be used to specify the character encoding of the source files
|
||||
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
|
||||
|
|
@ -896,7 +901,8 @@ EXCLUDE_SYMBOLS =
|
|||
# that contain example code fragments that are included (see the \include
|
||||
# command).
|
||||
|
||||
EXAMPLE_PATH =
|
||||
EXAMPLE_PATH = ../test/src \
|
||||
../test/dox
|
||||
|
||||
# If the value of the EXAMPLE_PATH tag contains directories, you can use the
|
||||
# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and
|
||||
|
|
@ -910,7 +916,7 @@ EXAMPLE_PATTERNS = *
|
|||
# irrespective of the value of the RECURSIVE tag.
|
||||
# The default value is: NO.
|
||||
|
||||
EXAMPLE_RECURSIVE = NO
|
||||
EXAMPLE_RECURSIVE = YES
|
||||
|
||||
# The IMAGE_PATH tag can be used to specify one or more files or directories
|
||||
# that contain images that are to be included in the documentation (see the
|
||||
|
|
@ -1582,7 +1588,7 @@ MATHJAX_CODEFILE =
|
|||
# The default value is: YES.
|
||||
# This tag requires that the tag GENERATE_HTML is set to YES.
|
||||
|
||||
SEARCHENGINE = YES
|
||||
SEARCHENGINE = NO
|
||||
|
||||
# When the SERVER_BASED_SEARCH tag is enabled the search engine will be
|
||||
# implemented using a web server instead of a web client using Javascript. There
|
||||
|
|
@ -2197,7 +2203,7 @@ DIA_PATH =
|
|||
# and usage relations if the target is undocumented or is not a class.
|
||||
# The default value is: YES.
|
||||
|
||||
HIDE_UNDOC_RELATIONS = YES
|
||||
HIDE_UNDOC_RELATIONS = NO
|
||||
|
||||
# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is
|
||||
# available from the path. This tool is part of Graphviz (see:
|
||||
|
|
@ -2206,7 +2212,7 @@ HIDE_UNDOC_RELATIONS = YES
|
|||
# set to NO
|
||||
# The default value is: YES.
|
||||
|
||||
HAVE_DOT = YES
|
||||
HAVE_DOT = NO
|
||||
|
||||
# The DOT_NUM_THREADS specifies the number of dot invocations doxygen is allowed
|
||||
# to run in parallel. When set to 0 doxygen will base this on the number of
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@
|
|||
\defgroup DFT_GROUP Алгоритмы дискретного и быстрого преобразования Фурье
|
||||
\ingroup SPECTRAL_GROUP
|
||||
|
||||
Алгоритмы дискретного и быстрого преобразования Фурье используют функции из библиотеки FFTW
|
||||
Алгоритмы дискретного и быстрого преобразования Фурье.
|
||||
|
||||
|
||||
\defgroup WIN_GROUP Функции оконного взвешивания
|
||||
|
|
|
|||
|
|
@ -7,7 +7,7 @@
|
|||
Функция рассчитывает \f$ n \f$-точечное дискретное преобразование Фурье
|
||||
вещественного сигнала \f$ x(m) \f$, \f$ m = 0 \ldots n-1 \f$.<BR>
|
||||
\f[
|
||||
Y(k) = \sum_{m = 0}^{n-1} x(m) \exp \left( -j \frac{2\pi}{n} m k \right),
|
||||
Y(k) = \sum_{m = 0}^{n-1} x(m) \cdot \exp \left( -j \cdot \frac{2\pi}{n} \cdot m \cdot k \right),
|
||||
\f]
|
||||
где \f$ k = 0 \ldots n-1 \f$.
|
||||
|
||||
|
|
@ -47,7 +47,7 @@
|
|||
Функция рассчитывает \f$ n \f$-точечное дискретное преобразование Фурье
|
||||
комплексного сигнала \f$ x(m) \f$, \f$ m = 0 \ldots n-1 \f$.<BR>
|
||||
\f[
|
||||
Y(k) = \sum_{m = 0}^{n-1} x(m) \exp \left( -j \frac{2\pi}{n} m k \right),
|
||||
Y(k) = \sum_{m = 0}^{n-1} x(m) \cdot \exp \left( -j \cdot \frac{2\pi}{n} \cdot m \cdot k \right),
|
||||
\f]
|
||||
где \f$ k = 0 \ldots n-1 \f$.
|
||||
|
||||
|
|
|
|||
|
|
@ -1,16 +1,133 @@
|
|||
/**************************************************************************************************
|
||||
Linspace array filling
|
||||
int linspace(double x0, double x1, int n, int type, double* x)
|
||||
***************************************************************************************************/
|
||||
/*! *************************************************************************************************
|
||||
\ingroup SPEC_MATH_COMMON_GROUP
|
||||
\fn int linspace(double x0, double x1, int n, int type, double* x)
|
||||
\brief Функция заполняет массив линейно-нарастающими, равноотстоящими значениями
|
||||
от `x0` до `x1`
|
||||
|
||||
Заполняет массив `x` длиной `n` значениями в диапазоне от \f$x_0\f$ до \f$x_1\f$.
|
||||
Функция поддерживает два типа заполнения в соответствии с параметром `type`:<BR>
|
||||
|
||||
Симметричное заполнение согласно выражению (параметр `type=DSPL_SYMMETRIC`):<BR>
|
||||
|
||||
\f$x(k) = x_0 + k \cdot dx\f$, here \f$dx = \frac{x_1 - x_0}{n-1}\f$, \f$k = 0 \ldots n-1.\f$
|
||||
|
||||
Периодическое заполнение (параметр `type=DSPL_PERIODIC`) согласно выражению:<BR>
|
||||
|
||||
\f$x(k) = x_0 + k \cdot dx\f$, here \f$dx = \frac{x_1 - x_0}{n}\f$, \f$k = 0 \ldots n-1.\f$
|
||||
|
||||
\param[in] x0 Начальное показателя \f$x_0\f$.<BR><BR>
|
||||
|
||||
\param[in] x1 Конечное значение \f$x_1\f$.<BR><BR>
|
||||
|
||||
\param[in] n Количество точек массива `x`.<BR><BR>
|
||||
|
||||
\param[in] type Тип заполнения:<BR>
|
||||
`DSPL_SYMMETRIC` - симметричное заполнение,<BR>
|
||||
`DSPL_PERIODIC` - периодическое заполнение.<BR><BR>
|
||||
|
||||
\param[in,out] x Указатель на вектор равноотстоящих значений .<BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Память должна быть выделена.<BR><BR>
|
||||
|
||||
\return
|
||||
`RES_OK` - функция выполнена успешно. <BR>
|
||||
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
|
||||
|
||||
\note
|
||||
Отличие периодического и симметричного заполнения можно понять из следующих примеров. <BR>
|
||||
Пример 1. Периодическое заполнение.
|
||||
\code
|
||||
double x[5];
|
||||
dspl_linspace(0, 5, 5, DSPL_PERIODIC, x);
|
||||
\endcode
|
||||
В массиве `x` будут лежать значения:
|
||||
\code
|
||||
0, 1, 2, 3, 4
|
||||
\endcode
|
||||
<BR><BR>
|
||||
Пример 2. Симметричное заполнение.
|
||||
\code
|
||||
double x[5];
|
||||
dspl_linspace(0, 5, 5, DSPL_SYMMETRIC, x);
|
||||
\endcode
|
||||
В массиве `x` будут лежать значения:
|
||||
\code
|
||||
0, 1.25, 2.5, 3.75, 5
|
||||
\endcode
|
||||
|
||||
\author
|
||||
Бахурин Сергей.
|
||||
www.dsplib.org
|
||||
|
||||
**************************************************************************************************** */
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
Logspace array filling
|
||||
int logspace(double x0, double x1, int n, int type, double* x)
|
||||
***************************************************************************************************/
|
||||
|
||||
|
||||
|
||||
/*! *************************************************************************************************
|
||||
\ingroup SPEC_MATH_COMMON_GROUP
|
||||
\fn int logspace(double x0, double x1, int n, int type, double* x)
|
||||
\brief Функция заполняет массив значениями логарифмической шкале
|
||||
|
||||
Заполняет массив `x` длиной `n` значениями в диапазоне от \f$10^{x_0}\f$ до \f$10^{x_1}\f$.<BR>
|
||||
Функция поддерживает два типа заполнения в соответствии с параметром `type`:<BR>
|
||||
|
||||
Симметричное заполнение согласно выражению:<BR>
|
||||
|
||||
\f$x(k) = 10^{x_0} \cdot dx^k\f$, где \f$dx = \sqrt[n-1]{10^{x_1 - x_0}}\f$,
|
||||
\f$k = 0 \ldots n-1.\f$
|
||||
|
||||
Периодическое заполнение согласно выражению:
|
||||
|
||||
\f$x(k) = 10^{x_0} \cdot dx^k\f$, где \f$dx = \sqrt[n]{10^{x_1 - x_0}}\f$,
|
||||
\f$k = 0 \ldots n-1.\f$<BR>
|
||||
|
||||
\param[in] x0 Начальное значение показателя \f$x_0\f$.<BR><BR>
|
||||
|
||||
\param[in] x1 Конечное значение показателя \f$x_1\f$.<BR><BR>
|
||||
|
||||
\param[in] n Количество точек массива `x`.<BR><BR>
|
||||
|
||||
\param[in] type Тип заполнения:<BR>
|
||||
`DSPL_SYMMETRIC` - симметричное заполнение,<BR>
|
||||
`DSPL_PERIODIC` - периодическое заполнение.<BR><BR>
|
||||
|
||||
\param[in,out] x Указатель на вектор значений в логарифмической шкале.<BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Память должна быть выделена.<BR><BR>
|
||||
|
||||
\return
|
||||
`RES_OK` - функция выполнена успешно. <BR>
|
||||
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
|
||||
|
||||
\note
|
||||
Отличие периодического и симметричного заполнения можно понять из следующих примеров. <BR>
|
||||
Пример 1. Периодическое заполнение.
|
||||
\code
|
||||
double x[5];
|
||||
dspl_logspace(-2, 3, 5, DSPL_PERIODIC, x);
|
||||
\endcode
|
||||
В массиве `x` будут лежать значения:
|
||||
\code
|
||||
0.01, 0.1, 1, 10, 100
|
||||
\endcode
|
||||
<BR><BR>
|
||||
Пример 2. Симметричное заполнение.
|
||||
\code
|
||||
double x[5];
|
||||
dspl_logspace(-2, 3, 5, DSPL_SYMMETRIC, x);
|
||||
\endcode
|
||||
В массиве `x` будут лежать значения:
|
||||
\code
|
||||
0.01 0.178 3.162 56.234 1000
|
||||
\endcode
|
||||
|
||||
\author
|
||||
Бахурин Сергей.
|
||||
www.dsplib.org
|
||||
|
||||
**************************************************************************************************** */
|
||||
|
||||
|
|
|
|||
|
|
@ -1,18 +1,106 @@
|
|||
|
||||
/*! *************************************************************************************************
|
||||
\ingroup SPEC_MATH_COMMON_GROUP
|
||||
\fn int polyval(double* a, int ord, double* x, int n, double* y)
|
||||
\brief Расчет вещественного полинома
|
||||
|
||||
Функция рассчитывает полином \f$P_N(x)\f$ \f$N-\f$ого порядка для вещественного
|
||||
аргумента, заданного вектором `x`.<BR>
|
||||
\f[
|
||||
P_N(x) = a_0 + a_1 \cdot x + a_2 \cdot x^2 + a_3 \cdot x^3 + ... a_N \cdot x^N.
|
||||
\f]
|
||||
Для расчета используется формула Горнера:<BR>
|
||||
\f[
|
||||
P_N(x) = a_0 + x \cdot (a_1 + x \cdot (a_2 + \cdot ( \ldots x \cdot (a_{N-1} + x\cdot a_N) \ldots )))
|
||||
\f]
|
||||
|
||||
\param[in] a Указатель на вектор вещественных коэффициентов полинома. <BR>
|
||||
Размер вектора `[ord+1 x 1]`.<BR>
|
||||
Коэффициент `a[0]` соответствует коэффициенту полинома \f$a_0\f$.<BR><BR>
|
||||
|
||||
\param[in] ord Порядок полинома \f$N\f$. <BR><BR>
|
||||
|
||||
\param[in] x Указатель на вектор аргумента полинома. <BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Значения полинома будут расчитаны для всех
|
||||
значений аргумента вектора `x`.<BR><BR>
|
||||
|
||||
\param[in] n Размер вектора агрумента полинома. <BR><BR>
|
||||
|
||||
\param[out] y Указатель на значения полинома для аргумента `x`. <BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Память должна быть выделена.<BR><BR>
|
||||
|
||||
/**************************************************************************************************
|
||||
Real polynomial evaluation
|
||||
\return
|
||||
`RES_OK` Полином расчитан успешно. <BR>
|
||||
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
|
||||
|
||||
\author
|
||||
Бахурин Сергей.
|
||||
www.dsplib.org
|
||||
|
||||
|
||||
int polyval(double* a, int ord, double* x, int n, double* y)
|
||||
***************************************************************************************************/
|
||||
**************************************************************************************************** */
|
||||
|
||||
|
||||
|
||||
|
||||
/*! *************************************************************************************************
|
||||
\ingroup SPEC_MATH_COMMON_GROUP
|
||||
\fn int polyval_cmplx(complex_t* a, int ord, complex_t* x, int n, complex_t* y)
|
||||
\brief Расчет комплексного полинома
|
||||
|
||||
Функция рассчитывает полином \f$P_N(x)\f$ \f$N-\f$ого порядка комплексного аргумента
|
||||
аргумента, заданного вектором `x`.<BR>
|
||||
\f[
|
||||
P_N(x) = a_0 + a_1 \cdot x + a_2 \cdot x^2 + a_3 \cdot x^3 + ... a_N \cdot x^N.
|
||||
\f]
|
||||
Для расчета используется формула Горнера:<BR>
|
||||
\f[
|
||||
P_N(x) = a_0 + x \cdot (a_1 + x \cdot (a_2 + \cdot ( \ldots x \cdot (a_{N-1} + x\cdot a_N) \ldots )))
|
||||
\f]
|
||||
|
||||
\param[in] aR Указатель на вектор реальной части коэффициентов полинома. <BR>
|
||||
Размер вектора `[ord+1 x 1]`.<BR>
|
||||
Коэффициент `aR[0]` соответствует коэффициенту полинома \f$a_0\f$.<BR><BR>
|
||||
|
||||
\param[in] aI Указатель на вектор мнимой части коэффициентов полинома. <BR>
|
||||
Размер вектора `[ord+1 x 1]`.<BR>
|
||||
Коэффициент `aI[0]` соответствует коэффициенту полинома \f$a_0\f$.<BR><BR>
|
||||
|
||||
\param[in] ord Порядок полинома \f$N\f$. <BR><BR>
|
||||
|
||||
\param[in] xR Указатель на вектор реальной части аргумента полинома. <BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Значения полинома будут расчитаны для всех
|
||||
значений аргумента вектора `x`.<BR><BR>
|
||||
|
||||
\param[in] xI Указатель на вектор мнимой части аргумента полинома. <BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Значения полинома будут расчитаны для всех
|
||||
значений аргумента вектора `x`.<BR><BR>
|
||||
|
||||
\param[in] n Размер вектора агрумента полинома. <BR><BR>
|
||||
|
||||
\param[out] yR Указатель вектор реальной части значения
|
||||
полинома для аргумента `x`. <BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Память должна быть выделена.<BR><BR>
|
||||
|
||||
\param[out] yI Указатель вектор реальной части значения
|
||||
полинома для аргумента `x`. <BR>
|
||||
Размер вектора `[n x 1]`.<BR>
|
||||
Память должна быть выделена.<BR><BR>
|
||||
|
||||
\return
|
||||
`RES_OK` Полином расчитан успешно. <BR>
|
||||
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
|
||||
|
||||
\author
|
||||
Бахурин Сергей.
|
||||
www.dsplib.org
|
||||
|
||||
**************************************************************************************************** */
|
||||
|
||||
/**************************************************************************************************
|
||||
Complex polynomial evaluation
|
||||
|
||||
int polyval_cmplx(complex_t* a, int ord, complex_t* x, int n, complex_t* y)
|
||||
***************************************************************************************************/
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,95 @@
|
|||
/*
|
||||
* Copyright (c) 2015-2018 Sergey Bakhurin
|
||||
* Digital Signal Processing Library [http://dsplib.org]
|
||||
*
|
||||
* This file is part of libdspl-2.0.
|
||||
*
|
||||
* is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU Lesser General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* DSPL is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License
|
||||
* along with Foobar. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "dspl.h"
|
||||
|
||||
int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2);
|
||||
int fft_p2(int n);
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
FFT bit reverse
|
||||
**************************************************************************************************/
|
||||
int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2)
|
||||
{
|
||||
static unsigned char rb_table[256] =
|
||||
{
|
||||
0x00, 0x80, 0x40, 0xC0, 0x20, 0xA0, 0x60, 0xE0, 0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,
|
||||
0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8, 0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,
|
||||
0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4, 0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,
|
||||
0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC, 0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,
|
||||
0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2, 0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,
|
||||
0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA, 0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,
|
||||
0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6, 0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,
|
||||
0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE, 0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,
|
||||
0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1, 0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,
|
||||
0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9, 0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,
|
||||
0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5, 0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,
|
||||
0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED, 0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,
|
||||
0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3, 0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,
|
||||
0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB, 0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,
|
||||
0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7, 0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,
|
||||
0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF, 0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF
|
||||
};
|
||||
|
||||
if(!x || !y)
|
||||
return ERROR_PTR;
|
||||
if(n<1 || p2 < 1)
|
||||
return ERROR_SIZE;
|
||||
|
||||
unsigned int v, c;
|
||||
|
||||
for(v = 0; v < n; v++)
|
||||
{
|
||||
c =
|
||||
((rb_table[ v & 0xff] << 24) |
|
||||
(rb_table[(v >> 8) & 0xff] << 16) |
|
||||
(rb_table[(v >> 16) & 0xff] << 8) |
|
||||
(rb_table[(v >> 24) & 0xff])) >> (32 - p2);
|
||||
RE(y[c]) = RE(x[v]);
|
||||
IM(y[c]) = IM(x[v]);
|
||||
}
|
||||
return RES_OK;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
FFT power of 2
|
||||
**************************************************************************************************/
|
||||
int fft_p2(int n)
|
||||
{
|
||||
int p = (n>0) ? n : 0;
|
||||
int p2 = 0;
|
||||
while((p = p>>1))
|
||||
p2++;
|
||||
if((1<<p2)!=n)
|
||||
return -1;
|
||||
return p2;
|
||||
}
|
||||
|
|
@ -0,0 +1,109 @@
|
|||
/*
|
||||
* Copyright (c) 2015-2018 Sergey Bakhurin
|
||||
* Digital Signal Processing Library [http://dsplib.org]
|
||||
*
|
||||
* This file is part of libdspl-2.0.
|
||||
*
|
||||
* is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU Lesser General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* DSPL is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License
|
||||
* along with Foobar. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "dspl.h"
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
Goertzel algorithm for real vector
|
||||
**************************************************************************************************/
|
||||
int DSPL_API goertzel(double *x, int n, int *ind, int k, complex_t *y)
|
||||
{
|
||||
|
||||
int m, p;
|
||||
double wR, wI;
|
||||
double alpha;
|
||||
double v[3];
|
||||
|
||||
if(!x || !y || !ind)
|
||||
return ERROR_PTR;
|
||||
|
||||
if(n < 1 || k < 1)
|
||||
return ERROR_SIZE;
|
||||
|
||||
for(p = 0; p < k; p++)
|
||||
{
|
||||
wR = cos(M_2PI * (double)ind[p] / (double)n);
|
||||
wI = sin(M_2PI * (double)ind[p] / (double)n);
|
||||
|
||||
alpha = 2.0 * wR;
|
||||
v[0] = v[1] = v[2] = 0.0;
|
||||
|
||||
for(m = 0; m < n; m++)
|
||||
{
|
||||
v[2] = v[1];
|
||||
v[1] = v[0];
|
||||
v[0] = x[m]+alpha*v[1] - v[2];
|
||||
}
|
||||
RE(y[p]) = wR * v[0] - v[1];
|
||||
IM(y[p]) = wI * v[0];
|
||||
}
|
||||
|
||||
return RES_OK;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
Goertzel algorithm for complex vector
|
||||
**************************************************************************************************/
|
||||
int DSPL_API goertzel_cmplx(complex_t *x, int n, int *ind, int k, complex_t *y)
|
||||
{
|
||||
|
||||
int m, p;
|
||||
complex_t w;
|
||||
double alpha;
|
||||
complex_t v[3];
|
||||
|
||||
if(!x || !y || !ind)
|
||||
return ERROR_PTR;
|
||||
|
||||
if(n < 1 || k < 1)
|
||||
return ERROR_SIZE;
|
||||
|
||||
for(p = 0; p < k; p++)
|
||||
{
|
||||
RE(w) = cos(M_2PI * (double)ind[p] / (double)n);
|
||||
IM(w) = sin(M_2PI * (double)ind[p] / (double)n);
|
||||
|
||||
alpha = 2.0 * RE(w);
|
||||
memset(v, 0, 3*sizeof(complex_t));
|
||||
|
||||
for(m = 0; m < n; m++)
|
||||
{
|
||||
RE(v[2]) = RE(v[1]);
|
||||
RE(v[1]) = RE(v[0]);
|
||||
RE(v[0]) = RE(x[m]) + alpha * RE(v[1]) - RE(v[2]);
|
||||
|
||||
IM(v[2]) = IM(v[1]);
|
||||
IM(v[1]) = IM(v[0]);
|
||||
IM(v[0]) = IM(x[m]) + alpha * IM(v[1]) - IM(v[2]);
|
||||
}
|
||||
|
||||
RE(y[p]) = CMRE(w, v[0]) - RE(v[1]);
|
||||
IM(y[p]) = CMIM(w, v[0]) - IM(v[1]);
|
||||
}
|
||||
|
||||
return RES_OK;
|
||||
}
|
||||
|
|
@ -0,0 +1,137 @@
|
|||
/*
|
||||
* Copyright (c) 2015-2018 Sergey Bakhurin
|
||||
* Digital Signal Processing Library [http://dsplib.org]
|
||||
*
|
||||
* This file is part of libdspl-2.0.
|
||||
*
|
||||
* is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU Lesser General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* DSPL is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License
|
||||
* along with Foobar. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "dspl.h"
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
Real polynomial evaluation
|
||||
***************************************************************************************************/
|
||||
int DSPL_API polyval(double* a, int ord, double* x, int n, double* y)
|
||||
{
|
||||
int k, m;
|
||||
|
||||
if(!a || !x || !y)
|
||||
return ERROR_PTR;
|
||||
if(ord<0)
|
||||
return ERROR_POLY_ORD;
|
||||
if(n<1)
|
||||
return ERROR_SIZE;
|
||||
|
||||
for(k = 0; k < n; k++)
|
||||
{
|
||||
y[k] = a[ord];
|
||||
for(m = ord-1; m>-1; m--)
|
||||
y[k] = y[k]*x[k] + a[m];
|
||||
}
|
||||
return RES_OK;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
Complex polynomial evaluation
|
||||
***************************************************************************************************/
|
||||
int DSPL_API polyval_cmplx(complex_t* a, int ord, complex_t* x, int n, complex_t* y)
|
||||
{
|
||||
int k, m;
|
||||
complex_t t;
|
||||
|
||||
if(!a || !x || !y)
|
||||
return ERROR_PTR;
|
||||
if(ord<0)
|
||||
return ERROR_POLY_ORD;
|
||||
if(n<1)
|
||||
return ERROR_SIZE;
|
||||
|
||||
for(k = 0; k < n; k++)
|
||||
{
|
||||
RE(y[k]) = RE(a[ord]);
|
||||
IM(y[k]) = IM(a[ord]);
|
||||
for(m = ord-1; m>-1; m--)
|
||||
{
|
||||
RE(t) = CMRE(y[k], x[k]);
|
||||
IM(t) = CMIM(y[k], x[k]);
|
||||
RE(y[k]) = RE(t) + RE(a[m]);
|
||||
IM(y[k]) = IM(t) + IM(a[m]);
|
||||
}
|
||||
}
|
||||
return RES_OK;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************************************************************************
|
||||
polynomial zeros to coeff reecalc
|
||||
***************************************************************************************************/
|
||||
int poly_z2a(complex_t *z, int nz, int ord, complex_t *a)
|
||||
{
|
||||
int k, ind, res;
|
||||
complex_t x[2];
|
||||
|
||||
if(!z || !a)
|
||||
return ERROR_PTR;
|
||||
if(ord<0)
|
||||
return ERROR_POLY_ORD;
|
||||
if(nz<1 || nz > ord)
|
||||
return ERROR_SIZE;
|
||||
|
||||
memset(a, 0, (ord+1) * sizeof(complex_t));
|
||||
RE(a[0]) = 1.0;
|
||||
IM(a[0]) = 0.0;
|
||||
|
||||
RE(x[1]) = 1.0;
|
||||
IM(x[1]) = 0.0;
|
||||
|
||||
|
||||
ind = 1;
|
||||
for(k=0; k < nz; k++)
|
||||
{
|
||||
RE(x[0]) = -RE(z[k]);
|
||||
IM(x[0]) = -IM(z[k]);
|
||||
res = conv_cmplx(a, ind, x, 2, a);
|
||||
if(res!=RES_OK)
|
||||
return res;
|
||||
ind++;
|
||||
}
|
||||
return RES_OK;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -42,8 +42,14 @@ p_conv_cmplx conv_cmplx ;
|
|||
p_dft dft ;
|
||||
p_dft_cmplx dft_cmplx ;
|
||||
p_filter_iir filter_iir ;
|
||||
p_goertzel goertzel ;
|
||||
p_goertzel_cmplx goertzel_cmplx ;
|
||||
p_linspace linspace ;
|
||||
p_logspace logspace ;
|
||||
p_polyval polyval ;
|
||||
p_polyval_cmplx polyval_cmplx ;
|
||||
|
||||
|
||||
#endif //BUILD_LIB
|
||||
|
||||
|
||||
|
|
@ -103,8 +109,13 @@ void* dspl_load()
|
|||
LOAD_FUNC(dft);
|
||||
LOAD_FUNC(dft_cmplx);
|
||||
LOAD_FUNC(filter_iir);
|
||||
LOAD_FUNC(goertzel);
|
||||
LOAD_FUNC(goertzel_cmplx);
|
||||
LOAD_FUNC(linspace);
|
||||
LOAD_FUNC(logspace);
|
||||
LOAD_FUNC(polyval);
|
||||
LOAD_FUNC(polyval_cmplx);
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -43,13 +43,16 @@ typedef double complex_t[2];
|
|||
|
||||
typedef struct
|
||||
{
|
||||
void *pfftw;
|
||||
complex_t *in;
|
||||
complex_t *out;
|
||||
size_t size;
|
||||
complex_t* w;
|
||||
complex_t* t0;
|
||||
complex_t* t1;
|
||||
int n;
|
||||
int p2;
|
||||
} fft_t;
|
||||
|
||||
|
||||
|
||||
|
||||
#define RE(x) (x[0])
|
||||
#define IM(x) (x[1])
|
||||
|
||||
|
|
@ -66,8 +69,6 @@ typedef struct
|
|||
|
||||
|
||||
|
||||
|
||||
|
||||
#define RES_OK 0
|
||||
|
||||
/* Error codes */
|
||||
|
|
@ -179,7 +180,7 @@ extern "C" {
|
|||
|
||||
#ifdef BUILD_LIB
|
||||
#define DECLARE_FUNC(type, fn, param)\
|
||||
type DSPL_API fn(param);\
|
||||
type DSPL_API fn(param);
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -198,8 +199,15 @@ DECLARE_FUNC(int, conv_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(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*);
|
||||
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*);
|
||||
DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*);
|
||||
DECLARE_FUNC(int, logspace, double COMMA double COMMA int COMMA int COMMA double*);
|
||||
DECLARE_FUNC(int, polyval, double* COMMA int COMMA double* COMMA int COMMA double*);
|
||||
DECLARE_FUNC(int, polyval_cmplx, complex_t* COMMA int COMMA complex_t* COMMA int COMMA complex_t*);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
|
|
|
|||
|
|
@ -42,8 +42,14 @@ p_conv_cmplx conv_cmplx ;
|
|||
p_dft dft ;
|
||||
p_dft_cmplx dft_cmplx ;
|
||||
p_filter_iir filter_iir ;
|
||||
p_goertzel goertzel ;
|
||||
p_goertzel_cmplx goertzel_cmplx ;
|
||||
p_linspace linspace ;
|
||||
p_logspace logspace ;
|
||||
p_polyval polyval ;
|
||||
p_polyval_cmplx polyval_cmplx ;
|
||||
|
||||
|
||||
#endif //BUILD_LIB
|
||||
|
||||
|
||||
|
|
@ -103,8 +109,13 @@ void* dspl_load()
|
|||
LOAD_FUNC(dft);
|
||||
LOAD_FUNC(dft_cmplx);
|
||||
LOAD_FUNC(filter_iir);
|
||||
LOAD_FUNC(goertzel);
|
||||
LOAD_FUNC(goertzel_cmplx);
|
||||
LOAD_FUNC(linspace);
|
||||
LOAD_FUNC(logspace);
|
||||
LOAD_FUNC(polyval);
|
||||
LOAD_FUNC(polyval_cmplx);
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -43,13 +43,16 @@ typedef double complex_t[2];
|
|||
|
||||
typedef struct
|
||||
{
|
||||
void *pfftw;
|
||||
complex_t *in;
|
||||
complex_t *out;
|
||||
size_t size;
|
||||
complex_t* w;
|
||||
complex_t* t0;
|
||||
complex_t* t1;
|
||||
int n;
|
||||
int p2;
|
||||
} fft_t;
|
||||
|
||||
|
||||
|
||||
|
||||
#define RE(x) (x[0])
|
||||
#define IM(x) (x[1])
|
||||
|
||||
|
|
@ -66,8 +69,6 @@ typedef struct
|
|||
|
||||
|
||||
|
||||
|
||||
|
||||
#define RES_OK 0
|
||||
|
||||
/* Error codes */
|
||||
|
|
@ -179,7 +180,7 @@ extern "C" {
|
|||
|
||||
#ifdef BUILD_LIB
|
||||
#define DECLARE_FUNC(type, fn, param)\
|
||||
type DSPL_API fn(param);\
|
||||
type DSPL_API fn(param);
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -198,8 +199,15 @@ DECLARE_FUNC(int, conv_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(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*);
|
||||
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*);
|
||||
DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*);
|
||||
DECLARE_FUNC(int, logspace, double COMMA double COMMA int COMMA int COMMA double*);
|
||||
DECLARE_FUNC(int, polyval, double* COMMA int COMMA double* COMMA int COMMA double*);
|
||||
DECLARE_FUNC(int, polyval_cmplx, complex_t* COMMA int COMMA complex_t* COMMA int COMMA complex_t*);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
|
|
|
|||
Plik binarny nie jest wyświetlany.
Ładowanie…
Reference in New Issue