mirror
/
libdspl-2.0
kopia lustrzana https://github.com/Dsplib/libdspl-2.0
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

1. Added function window 

2. Deleted test folder
3. Addeddded examples folder
4. Updated doxygen project
pull/6/head
Dsplib 2019-02-17 16:59:42 +03:00
rodzic b0b23ce5a7
commit 966d85edb4
68 zmienionych plików z 267 dodań i 2057 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

8
Makefile
Wyświetl plik

@ -22,10 +22,10 @@ all:
$(MAKE) -f Makefile.dspl
cp -r include/dspl.h release/include/dspl.h
cp -r include/dspl.c release/include/dspl.c
$(MAKE) -f Makefile.test
cp -r $(RELEASE_DIR)/$(DSPL_LIBNAME) test/bin/$(DSPL_LIBNAME)
$(MAKE) -f Makefile.verif
cp -r $(RELEASE_DIR)/$(DSPL_LIBNAME) verif/bin/$(DSPL_LIBNAME)
$(MAKE) -f Makefile.examples
cp -r $(RELEASE_DIR)/$(DSPL_LIBNAME) examples/bin/$(DSPL_LIBNAME)
@ -33,5 +33,7 @@ all:
clean:
$(MAKE) -f Makefile.dspl clean
$(MAKE) -f Makefile.test clean
$(MAKE) -f Makefile.verif clean
$(MAKE) -f Makefile.examples clean

56
Makefile.test → Makefile.examples
Wyświetl plik

@ -1,28 +1,28 @@
PRJ_DIR = test
SRC_DIR = $(PRJ_DIR)/src
BIN_DIR = $(PRJ_DIR)/bin
include Makefile.dirs
DSPL_C_FILE = $(INC_DIR)/dspl.c
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)
$(CC) $(CFLAGS) $< $(DSPL_O_FILE) -o $@ $(LFLAGS)
$(DSPL_O_FILE):$(DSPL_C_FILE)
$(CC) -c $(CFLAGS) $(DSPL_C_FILE) -o $(DSPL_O_FILE) $(LFLAGS)
clean:
rm -f $(BIN_DIR)/*.exe
rm -f $(BIN_DIR)/$(DSPL_LIBNAME)
rm -f $(DSPL_O_FILE)
PRJ_DIR = examples
SRC_DIR = $(PRJ_DIR)/src
BIN_DIR = $(PRJ_DIR)/bin
include Makefile.dirs
DSPL_C_FILE = $(INC_DIR)/dspl.c
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)
$(CC) $(CFLAGS) $< $(DSPL_O_FILE) -o $@ $(LFLAGS)
$(DSPL_O_FILE):$(DSPL_C_FILE)
$(CC) -c $(CFLAGS) $(DSPL_C_FILE) -o $(DSPL_O_FILE) $(LFLAGS)
clean:
rm -f $(BIN_DIR)/*.exe
rm -f $(BIN_DIR)/$(DSPL_LIBNAME)
rm -f $(DSPL_O_FILE)

20
dox/doxyfile_ru
Wyświetl plik

@ -32,7 +32,7 @@ DOXYFILE_ENCODING = UTF-8
# title of most generated pages and in a few other places.
# The default value is: My Project.
PROJECT_NAME = " libdspl-2.0"
PROJECT_NAME = libdspl-2.0
# The PROJECT_NUMBER tag can be used to enter a project or revision number. This
# could be handy for archiving the generated documentation or if some version
@ -44,7 +44,7 @@ PROJECT_NUMBER =
# for a project that appears at the top of each page and should give viewer a
# quick idea about the purpose of the project. Keep the description short.
PROJECT_BRIEF =
PROJECT_BRIEF = "Библиотека алгоритмов цифровой обработки сигналов"
# With the PROJECT_LOGO tag one can specify a logo or an icon that is included
# in the documentation. The maximum height of the logo should not exceed 55
@ -817,9 +817,9 @@ INPUT = ru \
../dspl/dox/ru \
../dspl/src \
../include \
../test/dox/ru \
../test/src \
../test/bin/gnuplot
../examples/src \
../examples/bin/gnuplot \
../examples/bin/img
# 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
@ -935,9 +935,9 @@ EXCLUDE_SYMBOLS =
# that contain example code fragments that are included (see the \include
# command).
EXAMPLE_PATH = ../test/src \
../test/dox \
../test/bin/gnuplot
EXAMPLE_PATH = ../examples/src \
../examples/bin/gnuplot \
../examples/bin/img
# 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
@ -958,8 +958,8 @@ EXAMPLE_RECURSIVE = YES
# that contain images that are to be included in the documentation (see the
# \image command).
IMAGE_PATH = ../test/bin/img \
ru/img
IMAGE_PATH = ru/img \
../examples/bin/img
# The INPUT_FILTER tag can be used to specify a program that doxygen should
# invoke to filter for each input file. Doxygen will invoke the filter program

11
dox/makedoc.sh
Wyświetl plik

@ -1,13 +1,14 @@
#!/bin/bash
find ../../ru -name "*.c" -exec cp -rf {} ../test/src \;
find ../../ru -name "*.dox" -exec cp -rf {} ../test/dox/ru \;
find ../../ru -name "*.plt" -exec cp -rf {} ../test/bin/gnuplot \;
#find ../../ru -name "*.c" -exec cp -rf {} ../test/src \;
#find ../../ru -name "*.dox" -exec cp -rf {} ../test/dox/ru \;
#find ../../ru -name "*.plt" -exec cp -rf {} ../test/bin/gnuplot \;
cd ../
mingw32-make
cd test/bin
cd examples/bin
for file in *.exe
do
"./$file"
@ -20,5 +21,5 @@ cd ../
mingw32-make clean
cd dox
pkill -x gnuplot
#pkill -x gnuplot

109
dspl/dox/ru/fourier_series.dox
Wyświetl plik

@ -1,48 +1,44 @@
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn int fourier_series_dec(double* t, double* s, int nt,
double period, int nw, double* w, complex_t* y)
double period, int nw, double* w, complex_t* y)
\brief Расчет коэффициентов разложения в ряд Фурье
Функция рассчитывает спектр периодического сигнала при усечении ряда Фурье<BR>
\param[in] t Указатель на массив моментов времени дискретизации
исходного сигнала `s`<BR>
размер вектора вектора `[nt x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] t Указатель на массив моментов времени дискретизации
исходного сигнала `s`<BR>
размер вектора вектора `[nt x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] s Указатель на массив значений исходного сигнала`s`<BR>
размер вектора вектора `[nt x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] s Указатель на массив значений исходного сигнала`s`.<BR>
Размер вектора `[nt x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] nt размер выборки исходного сигнала.<BR>
Значение должно быть положительным.<BR><BR>
\param[in] nt Размер выборки исходного сигнала.<BR>
Значение должно быть положительным.<BR><BR>
\param[in] period Период повторения сигнала.<BR><BR>
\param[in] period Период повторения сигнала.<BR><BR>
\param[in] nw Размер усеченного ряда Фурье.<BR><BR>
\param[in] nw Размер усеченного ряда Фурье.<BR><BR>
\param[out] w Указатель на массив частот спектра
периодического сигнала.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[out] w Указатель на массив частот спектра
периодического сигнала.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[out] y Указатель массив комплексных значений спектра
периодического сигнала.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[out] y Указатель массив комплексных значений спектра
периодического сигнала.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` Массивы нулей и полюсов рассчитаны успешно.<BR>
В противном случае
\ref ERROR_CODE_GROUP "код ошибки".<BR>
`RES_OK` Коэффициенты ряда Фурье рассчитаны успешно.<BR>
В противном случае
\ref ERROR_CODE_GROUP "код ошибки".<BR>
\note
Для расчета спектра сигнала используется численное интегрирование
@ -53,8 +49,8 @@
<BR>
\author
Бахурин Сергей
www.dsplib.org
Бахурин Сергей
www.dsplib.org
***************************************************************************** */
@ -66,7 +62,7 @@
/*! ****************************************************************************
\ingroup DFT_GROUP
\fn int fourier_series_rec(double* w, complex_t* s, int nw,
double *t, int nt, complex_t* y)
double *t, int nt, complex_t* y)
\brief Восстановление сигнала при усечении ряда Фурье
@ -76,41 +72,36 @@
s(t) = \sum\limits_{n = 0}^{n_{\omega}-1} S(\omega_n) \exp(j\omega_n t)
\f]
\param[in] w Указатель на массив частот \f$\omega_n\f$
усеченного ряда Фурье.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена и заполнена.<BR><BR>
\param[in] w Указатель на массив частот \f$\omega_n\f$
усеченного ряда Фурье.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена и заполнена.<BR><BR>
\param[in] s Указатель на массив значений спектра
\f$S(\omega_n)\f$.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена и заполнена.<BR><BR>
\param[in] s Указатель на массив значений спектра
\f$S(\omega_n)\f$.<BR>
Размер вектора `[nw x 1]`.<BR>
Память должна быть выделена и заполнена.<BR><BR>
\param[in] nw Количество членов усеченного ряда Фурье.<BR>
Значение должно быть положительным.<BR><BR>
\param[in] nw Количество членов усеченного ряда Фурье.<BR>
Значение должно быть положительным.<BR><BR>
\param[in] t Указатель на массив временных отсчетов
восстановленного сигнала.<BR>
Размер вектора `[nt x 1]`.<BR>
Память должна быть выделена и заполнена.<BR><BR>
<BR><BR>
\param[in] t Указатель на массив временных отсчетов
восстановленного сигнала.<BR>
Размер вектора `[nt x 1]`.<BR>
Память должна быть выделена и заполнена.<BR><BR>
<BR><BR>
\param[in] nt Размер вектора времени и восстановленного сигнала.
<BR><BR>
\param[in] nt Размер вектора времени и восстановленного сигнала.<BR><BR>
\param[out] y Указатель на массив восстановленного сигнала.<BR>
Размер вектора `[nt x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[out] y Указатель на массив восстановленного сигнала.<BR>
Размер вектора `[nt x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` Массивы нулей и полюсов рассчитаны успешно.<BR>
В противном случае
\ref ERROR_CODE_GROUP "код ошибки".<BR>
`RES_OK` Массивы нулей и полюсов рассчитаны успешно.<BR>
В противном случае
\ref ERROR_CODE_GROUP "код ошибки".<BR>
\note
Выходной восстановленный сигнал в общем случае является комплексным.
@ -121,7 +112,7 @@ s(t) = \sum\limits_{n = 0}^{n_{\omega}-1} S(\omega_n) \exp(j\omega_n t)
<BR>
\author
Бахурин Сергей
www.dsplib.org
Бахурин Сергей
www.dsplib.org
***************************************************************************** */

66
dspl/dox/ru/win.dox 100644
Wyświetl plik

@ -0,0 +1,66 @@
/*! ****************************************************************************
\ingroup WIN_GROUP
\fn int window(double* w, int n, int win_type, double param)
\brief Расчет функции оконного взвешивания
Функция рассчитывает периодическую или симметричную оконную функцию
в соответствии с параметром `win_type`.<BR>
Периодическая оконная функция используется для спектрального анализа,
а симметричная оконная функция может быть использована для синтеза
КИХ-фильтров.<BR>
\param [in,out] w Указатель на вектор оконной функции.<BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR>
Рассчитанная оконная функция будет
помещена по данному адресу.<BR><BR>
\param [in] n Размер вектора `w` оконной функции.<BR><BR>
\param [in] win_type Комбинация флагов для задания типа оконной
функции.<BR>
Для задания типа окна используется
комбинация битовых масок
`DSPL_WIN_MASK | DSPL_WIN_SYM_MASK`.<BR>
Маска `DSPL_WIN_MASK` задает тип оконной функции.
Может принимать следующие значения:<BR>
Значение `DSPL_WIN_MASK` | Описание
----------------------------|----------------------------
`DSPL_WIN_BARTLETT` |Непараметрическое окно Бартлетта
`DSPL_WIN_BARTLETT_HANN` |Непараметрическое окно Бартлетта-Ханна
`DSPL_WIN_BLACKMAN` |Непараметрическое окно Блэкмана
`DSPL_WIN_BLACKMAN_HARRIS` |Непараметрическое окно Блэкмана-Харриса
`DSPL_WIN_BLACKMAN_NUTTALL` |Непараметрическое окно Блэкмана-Натталла
`DSPL_WIN_CHEBY` |Параметрическое окно Дольф-Чебышева. Данное окно всегда является симметричным и игнорирует параметр `DSPL_WIN_SYM_MASK`. Параметр `param` задает уровень боковых лепестков в дБ.
`DSPL_WIN_COS` |Непараметрическое косинус-окно
`DSPL_WIN_FLAT_TOP` |Непараметрическое окно с максимально плоской вершиной
`DSPL_WIN_GAUSSIAN` |Параметрическое окно Гаусса
`DSPL_WIN_HAMMING` |Непараметрическое окно Хемминга
`DSPL_WIN_HANN` |Непараметрическое окно Ханна
`DSPL_WIN_LANCZOS` |Непараметрическое окно Ланкзоса
`DSPL_WIN_NUTTALL` |Непараметрическое окно Натталла
`DSPL_WIN_RECT` |Непараметрическое прямоугольное окно
<BR><BR>Маска `DSPL_WIN_SYM_MASK` задает симметричное
или периодическое окно:<BR><BR>
Значение `DSPL_WIN_SYM_MASK` | Описание
--------------------------------|----------------------------
`DSPL_WIN_SYMMETRIC` |Симметричное окно (по умолчанию)
`DSPL_WIN_PERIODIC` |Периодическое окно
<BR><BR>
\param [in] param Параметр окна.<BR>
Данный параметр применяется только для параметрических оконных функций. <BR>
Для непараметрических окон игнорируется.<BR><BR>
\return
`DSPL_OK` если оконная функция рассчитана успешно.<BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
\author
Бахурин Сергей.
www.dsplib.org
***************************************************************************** */

25
dspl/src/dspl_internal.h
Wyświetl plik

@ -69,20 +69,21 @@ void dft16 (complex_t *x, complex_t* y);
void transpose4x4(complex_t *x, complex_t* y);
int win_bartlett (double *w, int n, int win_type);
/* Window functions */
int win_bartlett (double *w, int n, int win_type);
int win_bartlett_hann (double *w, int n, int win_type);
int win_blackman (double *w, int n, int win_type);
int win_blackman_harris (double *w, int n, int win_type);
int win_blackman (double *w, int n, int win_type);
int win_blackman_harris (double *w, int n, int win_type);
int win_blackman_nuttall(double *w, int n, int win_type);
int win_cos (double *w, int n, int win_type);
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_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);
int win_cheby (double *w, int n, double param);
int win_cos (double *w, int n, int win_type);
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_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);

67
dspl/src/win.c
Wyświetl plik

@ -45,6 +45,8 @@ int window(double* w, int n, int win_type, double param)
return win_blackman_harris(w, n, win_type);
case DSPL_WIN_BLACKMAN_NUTTALL:
return win_blackman_nuttall(w, n, win_type);
case DSPL_WIN_CHEBY:
return win_cheby(w, n, param);
case DSPL_WIN_FLAT_TOP:
return win_flat_top(w, n, win_type);
case DSPL_WIN_GAUSSIAN:
@ -60,9 +62,9 @@ int window(double* w, int n, int win_type, double param)
case DSPL_WIN_RECT:
return win_rect(w, n);
case DSPL_WIN_COS:
return win_cos(w, n, win_type);
return win_cos(w, n, win_type);
default:
return ERROR_WIN_TYPE;
return ERROR_WIN_TYPE;
}
return RES_OK;
}
@ -166,6 +168,67 @@ 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;
}
/******************************************************************************

0
test/bin/.gitignore → examples/bin/.gitignore vendored
Wyświetl plik

0
test/bin/dat/.gitignore → examples/bin/dat/.gitignore vendored
Wyświetl plik

0
test/bin/gnuplot/butter_ap_test.plt → examples/bin/gnuplot/butter_ap_test.plt
Wyświetl plik

0
test/bin/gnuplot/cheby1_ap_test.plt → examples/bin/gnuplot/cheby1_ap_test.plt
Wyświetl plik

0
test/bin/gnuplot/sinc_test.plt → examples/bin/gnuplot/sinc_test.plt
Wyświetl plik

0
test/bin/gnuplot/sine_int.plt → examples/bin/gnuplot/sine_int.plt
Wyświetl plik

0
test/bin/img/.gitignore → examples/bin/img/.gitignore vendored
Wyświetl plik

0
test/obj/.gitignore → examples/obj/.gitignore vendored
Wyświetl plik

0
test/src/butter_ap_test.c → examples/src/butter_ap_test.c
Wyświetl plik

0
test/src/cheby1_ap_test.c → examples/src/cheby1_ap_test.c
Wyświetl plik

0
test/src/cheby2_ap_zp_test.c → examples/src/cheby2_ap_zp_test.c
Wyświetl plik

0
test/src/conv_fft_cmplx_test.c → examples/src/conv_fft_cmplx_test.c
Wyświetl plik

0
test/src/dft_cmplx_test.c → examples/src/dft_cmplx_test.c
Wyświetl plik

0
test/src/dft_test.c → examples/src/dft_test.c
Wyświetl plik

0
test/src/dspl_info_test.c → examples/src/dspl_info_test.c
Wyświetl plik

0
test/src/ellip_landen_test.c → examples/src/ellip_landen_test.c
Wyświetl plik

0
test/src/fft_cmplx_test.c → examples/src/fft_cmplx_test.c
Wyświetl plik

0
test/src/fft_test.c → examples/src/fft_test.c
Wyświetl plik

0
test/src/matrix_print.c → examples/src/matrix_print.c
Wyświetl plik

0
test/src/matrix_transpose.c → examples/src/matrix_transpose.c
Wyświetl plik

0
test/src/sinc_test.c → examples/src/sinc_test.c
Wyświetl plik

0
test/src/sine_int_test.c → examples/src/sine_int_test.c
Wyświetl plik

2
include/dspl.c
Wyświetl plik

@ -125,6 +125,7 @@ p_trapint_cmplx trapint_cmplx ;
p_unwrap unwrap ;
p_verif verif ;
p_verif_cmplx verif_cmplx ;
p_window window ;
p_writebin writebin ;
p_writetxt writetxt ;
p_writetxt_3d writetxt_3d ;
@ -273,6 +274,7 @@ void* dspl_load()
LOAD_FUNC(unwrap);
LOAD_FUNC(verif);
LOAD_FUNC(verif_cmplx);
LOAD_FUNC(window);
LOAD_FUNC(writebin);
LOAD_FUNC(writetxt);
LOAD_FUNC(writetxt_3d);

13
include/dspl.h
Wyświetl plik

@ -109,7 +109,7 @@ typedef struct
/* K 0x11xxxxxx*/
/* L 0x12xxxxxx*/
/* M 0x13xxxxxx*/
#define ERROR_MATRIX_SIZE 0x13011909
#define ERROR_MATRIX_SIZE 0x13011909
/* N 0x14xxxxxx*/
#define ERROR_NEGATIVE 0x14050701
/* O 0x15xxxxxx*/
@ -129,8 +129,9 @@ typedef struct
#define ERROR_UNWRAP 0x21142318
/* V 0x22xxxxxx*/
/* W 0x23xxxxxx*/
#define ERROR_WIN_TYPE 0x23092025
#define ERROR_WIN_PARAM 0x23091601
#define ERROR_WIN_SYM 0x23091925
#define ERROR_WIN_TYPE 0x23092025
/* X 0x24xxxxxx*/
/* Y 0x25xxxxxx*/
/* Z 0x26xxxxxx*/
@ -150,7 +151,7 @@ typedef struct
#define DSPL_FLAG_UNWRAP 0x00000002
#define DSPL_WIN_SYM_MASK 0x00000001
#define DSPL_WIN_MASK 0x000FFFFE
#define DSPL_WIN_MASK 0x00FFFFFE
#define DSPL_WIN_SYMMETRIC DSPL_SYMMETRIC
#define DSPL_WIN_PERIODIC DSPL_PERIODIC
@ -169,6 +170,7 @@ typedef struct
#define DSPL_WIN_NUTTALL 0x00008000
#define DSPL_WIN_RECT 0x00010000
#define DSPL_WIN_COS 0x00040000
#define DSPL_WIN_CHEBY 0x00080000
#define ELLIP_ITER 16
@ -705,6 +707,11 @@ DECLARE_FUNC(int, verif_cmplx, complex_t* x
COMMA double eps
COMMA double* err);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, window, double* w
COMMA int n
COMMA int win_type
COMMA double param);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, writebin, void*
COMMA int
COMMA int

2
release/include/dspl.c
Wyświetl plik

@ -125,6 +125,7 @@ p_trapint_cmplx trapint_cmplx ;
p_unwrap unwrap ;
p_verif verif ;
p_verif_cmplx verif_cmplx ;
p_window window ;
p_writebin writebin ;
p_writetxt writetxt ;
p_writetxt_3d writetxt_3d ;
@ -273,6 +274,7 @@ void* dspl_load()
LOAD_FUNC(unwrap);
LOAD_FUNC(verif);
LOAD_FUNC(verif_cmplx);
LOAD_FUNC(window);
LOAD_FUNC(writebin);
LOAD_FUNC(writetxt);
LOAD_FUNC(writetxt_3d);

13
release/include/dspl.h
Wyświetl plik

@ -109,7 +109,7 @@ typedef struct
/* K 0x11xxxxxx*/
/* L 0x12xxxxxx*/
/* M 0x13xxxxxx*/
#define ERROR_MATRIX_SIZE 0x13011909
#define ERROR_MATRIX_SIZE 0x13011909
/* N 0x14xxxxxx*/
#define ERROR_NEGATIVE 0x14050701
/* O 0x15xxxxxx*/
@ -129,8 +129,9 @@ typedef struct
#define ERROR_UNWRAP 0x21142318
/* V 0x22xxxxxx*/
/* W 0x23xxxxxx*/
#define ERROR_WIN_TYPE 0x23092025
#define ERROR_WIN_PARAM 0x23091601
#define ERROR_WIN_SYM 0x23091925
#define ERROR_WIN_TYPE 0x23092025
/* X 0x24xxxxxx*/
/* Y 0x25xxxxxx*/
/* Z 0x26xxxxxx*/
@ -150,7 +151,7 @@ typedef struct
#define DSPL_FLAG_UNWRAP 0x00000002
#define DSPL_WIN_SYM_MASK 0x00000001
#define DSPL_WIN_MASK 0x000FFFFE
#define DSPL_WIN_MASK 0x00FFFFFE
#define DSPL_WIN_SYMMETRIC DSPL_SYMMETRIC
#define DSPL_WIN_PERIODIC DSPL_PERIODIC
@ -169,6 +170,7 @@ typedef struct
#define DSPL_WIN_NUTTALL 0x00008000
#define DSPL_WIN_RECT 0x00010000
#define DSPL_WIN_COS 0x00040000
#define DSPL_WIN_CHEBY 0x00080000
#define ELLIP_ITER 16
@ -705,6 +707,11 @@ DECLARE_FUNC(int, verif_cmplx, complex_t* x
COMMA double eps
COMMA double* err);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, window, double* w
COMMA int n
COMMA int win_type
COMMA double param);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, writebin, void*
COMMA int
COMMA int

21
test/bin/gnuplot/fourier_series_pimp.plt
Wyświetl plik

@ -1,21 +0,0 @@
set terminal pngcairo size 560,480 enhanced font 'Verdana,8'
set output 'img/fourier_series_pimp.png'
unset key
set grid
set lmargin 8
set multiplot layout 2,1 rowsfirst
set xlabel 'Частота, рад/с'
#
set ylabel 'Амплитудный спектр'
plot[-10*pi:10*pi] 'dat/fourier_series_pimp_mag.txt' with impulses lt 1 ,\
'dat/fourier_series_pimp_mag.txt' with points pt 7 ps 0.5 lt 1
#
set ylabel 'Фазовый спектр'
plot[-10*pi:10*pi] 'dat/fourier_series_pimp_phi.txt' with impulses lt 1 ,\
'dat/fourier_series_pimp_phi.txt' with points pt 7 ps 0.5 lt 1
unset multiplot

37
test/bin/gnuplot/fourier_series_rec.plt
Wyświetl plik

@ -1,37 +0,0 @@
### Start multiplot (2x2 layout)
set terminal pngcairo size 800,640 enhanced font 'Verdana,8'
set output 'img/fourier_series_rec.png'
unset key
set multiplot layout 4,2 rowsfirst
#
plot 'dat/fourier_series_pimp0.txt' with lines,\
'dat/fourier_series_pimp_rec_5.txt' with lines
#
plot 'dat/fourier_series_saw0.txt' with lines,\
'dat/fourier_series_saw_rec_5.txt' with lines
#
plot 'dat/fourier_series_pimp0.txt' with lines,\
'dat/fourier_series_pimp_rec_9.txt' with lines
#
plot 'dat/fourier_series_saw0.txt' with lines,\
'dat/fourier_series_saw_rec_9.txt' with lines
#
plot 'dat/fourier_series_pimp0.txt' with lines,\
'dat/fourier_series_pimp_rec_21.txt' with lines
#
plot 'dat/fourier_series_saw0.txt' with lines,\
'dat/fourier_series_saw_rec_21.txt' with lines
#
plot 'dat/fourier_series_pimp0.txt' with lines,\
'dat/fourier_series_pimp_rec_61.txt' with lines
#
plot 'dat/fourier_series_saw0.txt' with lines,\
'dat/fourier_series_saw_rec_61.txt' with lines
unset multiplot

24
test/bin/gnuplot/fourier_transform_ex_gauss.plt
Wyświetl plik

@ -1,24 +0,0 @@
set terminal pngcairo size 560,480 enhanced font 'Verdana,8'
set output 'img/fourier_series_pimp.png'
unset key
set grid
set lmargin 8
set multiplot layout 2,1 rowsfirst
#
set xlabel 't, c'
set ylabel 's(t)'
plot[-10:10] 'dat/fourier_transform_ex_gauss_time_0.5.txt' with lines,\
'dat/fourier_transform_ex_gauss_time_1.0.txt' with lines,\
'dat/fourier_transform_ex_gauss_time_2.0.txt' with lines
#
set xlabel 'w, рад/c'
set ylabel 'S(w)'
plot[-4*pi:4*pi] 'dat/fourier_transform_ex_gauss_freq_0.5.txt' with lines,\
'dat/fourier_transform_ex_gauss_freq_1.0.txt' with lines,\
'dat/fourier_transform_ex_gauss_freq_2.0.txt' with lines
unset multiplot

21
test/dox/ru/dft_indexation.dox
Wyświetl plik

@ -1,21 +0,0 @@
/*!
\example dft_indexation.c
\brief Расчет данных для построения рисунка 1 и рисунка 5 раздела
<a href = "http://ru.dsplib.org/content/dft_freq/dft_freq.html">
"Индексация и перестановка спектральных отсчетов дискретного преобразования Фурье"
</a>.
<BR>
Программа рассчитает \f$ N = 30 \f$ точечное ДПФ сигнала
\f[
s(n) = \exp \left( 2 \pi n \frac{f_0}{F_{\textrm{s}}} \right),~~ n = 0 \ldots N-1,
\f]
а также производит перестановку спектральных отсчетов полученного ДПФ при использовании
функции `fft_shift`.<BR>
Данные расчета сохраняются в текстовые файлы для построения графиков с
использованием пакетов tikz и pgfplots системы LaTeX.<BR>
Для построения графиков при помощи пакета GNUPlot Необходимо выполнить скрипт
*/

35
test/dox/ru/filter_approx.dox
Wyświetl plik

@ -1,35 +0,0 @@
/*!
\example filter_approx.c
\brief Расчет данных для построения графиков раздела
«<a href = "http://ru.dsplib.org/content/filter_approximation/filter_approximation.html">
Расчет аналогового фильтра. Постановка задачи и способы аппроксимации АЧХ нормированного ФНЧ
</a>».
Программа сохраняет текстовые файлы данных для построения графиков аппроксимирующих функций \f$ F_N^2(\omega)\f$ и
квадрата АЧХ \f$ |H(j\omega)|^2\f$ нормированных аналоговых фильтров нижних частот:
<dl>
<dt>`dat/butter_r.txt`</dt>
<dd>
Файл значений квадрата аппроксимирующей функции фильтра Баттерворта
</dd>
<dt>`dat/butter_h.txt`</dt>
<dd>
Файл значений квадрата АЧХ фильтра Баттерворта в линейном масштабе
</dd>
<dt>`dat/butter_hdb.txt`</dt>
<dd>
Файл значений квадрата АЧХ фильтра Баттерворта в логарифмическом масштабе (дБ)
</dd>
</dl>
*/

24
test/dox/ru/fourier_series_pimp_spectrum.dox
Wyświetl plik

@ -1,24 +0,0 @@
/*! ****************************************************************************
\example fourier_series_pimp_spectrum.c
Расчет амплитудного и фазового спектра периодической последовательности
прямоугольных импульсов
<BR>
Программа формирует один период последовательности прямоугольных импульсов
и производит расчет амплитудного и фазового спектра.
Результаты расчета сохраняются в файлы и выводятся на график программой GNUPlot
Скрипт GNUPLOT `fourier_series_pimp.plt`
для построения графиков из текстовых файлов:
\include fourier_series_pimp.plt
График будет отображен на экране и сохранен в файл `img/fourier_series_pimp.png`
\image html fourier_series_pimp.png
***************************************************************************** */

33
test/dox/ru/fourier_series_rec.dox
Wyświetl plik

@ -1,33 +0,0 @@
/*! ****************************************************************************
\example fourier_series_rec.c
Программа производит расчет данных для представления
периодической последовательности прямоугольных импульсов и пилообразного сигнала
усеченным рядом Фурье.
<BR>
Программа формирует `P=3` периодов последовательности прямоугольных импульсов
и пилообразного сигнала и производит их аппроксимацию усеченным рядом Фурье в
комплексной форме при различном количестве членов усеченного ряда: от 2 до 30.
Для каждого усеченного ряда производится расчет представления периодической
последовательности прямоугольных импульсов и пилообразного сигнала и сохранение
полученных данных в текстовые файлы для построения графиков.
Программа GNUPLOT произведет построение графика
по сохраненным в файлах данным.
Скрипт GNUPLOT `fourier_series_rec.plt`
для построения графиков из текстовых файлов:
\include fourier_series_rec.plt
График будет отображен на экране и сохранен в файл `img/fourier_series_rec.png`
\image html fourier_series_rec.png
***************************************************************************** */

28
test/dox/ru/sinc_test.dox
Wyświetl plik

@ -1,28 +0,0 @@
/*! ****************************************************************************
\example sinc_test.c Пример использования функции \ref sinc
Программа производит расчет функции
\f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$
для 3 различных значений параметра \f$ a\f$.
<BR>
В каталоге `dat` будут созданы три файла:<BR>
<pre>
sinc_test_1.0.txt параметр a = 1.0
sinc_test_pi.txt параметр a = pi
sinc_test_2pi.txt параметр a = 2*pi
</pre>
Кроме того программа GNUPLOT произведет построение графика
по сохраненным в файлах данным.
Скрипт GNUPLOT `sinc_test.plt` для построения графиков из текстовых файлов:
\include sinc_test.plt
График будет отображен на экране и сохранен в файл `img/sinc_test.png`
\image html sinc_test.png
***************************************************************************** */

31
test/octave/fft_n1n2.m
Wyświetl plik

@ -1,31 +0,0 @@
clear all;
close all;
clc;
N = 16;
N1 = 4;
N2 = 4;
% входной сигнал это вектор столбец размерности [N x 1]
x = (0:N-1)';
%x = [4;6;8;10];
for n1 = 0 : N1-1
for k2 = 0 : N2-1
W(n1 + 1, k2 + 1) = exp(-2i * pi * n1 * k2 / N);
end
end
A = reshape(x, N2, N1);
B = A.';
D = fft(B);
F = D.*W;
G = F.';
H = fft(G);
P = H.';
y = [reshape(P,N,1), fft(x)]

20
test/octave/fft_performance.m
Wyświetl plik

@ -1,20 +0,0 @@
clear all; close all; clc;
n = 8388608;
m = 4;
x0 = (0:n-1)+1i*(0:n-1);
while(n > 4)
x = x0(1:n);
fprintf('n = %12d ', n);
t0 = tic();
for k = 1:m
X = fft(x);
end
dt = toc(t0) / m;
fprintf('%12.6f\n', dt*1000);
n = n/2;
m = m*1.5;
end

22
test/octave/fft_twiddle.m
Wyświetl plik

@ -1,22 +0,0 @@
clear all; close all; clc;
u = - 2*pi/7;
w = [ ( cos(u) + cos(2*u) + cos(3*u)) / 3 - 1;
(2*cos(u) - cos(2*u) - cos(3*u)) / 3;
( cos(u) -2*cos(2*u) + cos(3*u)) / 3;
( cos(u) + cos(2*u) - 2*cos(3*u)) / 3;
( sin(u) + sin(2*u) - sin(3*u)) / 3;
(2*sin(u) - sin(2*u) + sin(3*u)) / 3;
( sin(u) -2*sin(2*u) - sin(3*u)) / 3;
( sin(u) + sin(2*u) + 2*sin(3*u)) / 3;]
fprintf('DFT 7 coeff\n');
for k = 1:length(w)
fprintf('#define DFT7_W%d % -.24f\n', k, w(k));
end

46
test/src/butter_ap.c
Wyświetl plik

@ -1,46 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define ORD 3
#define N 1000
#define K 1024
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1];
double Rp = 1.0;
double w[N], mag[N], phi[N], tau[N];
double t[K], h[K];
fft_t pfft;
int k;
int res = butter_ap(Rp, ORD, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD+1; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
freqs_resp(b, a, ORD, w, N, DSPL_FLAG_LOG|DSPL_FLAG_UNWRAP, mag, phi, tau);
writetxt(w, mag, N, "dat/butter_ap_test_mag.txt");
writetxt(w, phi, N, "dat/butter_ap_test_phi.txt");
writetxt(w, tau, N, "dat/butter_ap_test_tau.txt");
memset(&pfft, 0, sizeof(fft_t));
freqs2time(b, a, ORD, 200.0, K, &pfft, t,h);
writetxt(t, h, K, "dat/butter_ap_test_time.txt");
dspl_free(handle); // free dspl handle
return 0;
}

47
test/src/cheby1_ap_approx.c
Wyświetl plik

@ -1,47 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 1000
int main(int argc, char* argv[])
{
double w[N], r[N], h[N];
double ep2, Gp2;
int ord, k;
void* handle;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(0, 2.5, N, DSPL_PERIODIC, w);
ord = 3;
Gp2 = 0.9;
ep2 = 1.0 / Gp2 -1;
cheby_poly1(w, N, ord, r);
for(k = 0; k < N; k++)
{
r[k] *= r[k];
h[k] = 6.0/(1.0 + ep2*r[k]);
w[k] *= 4;
}
writetxt(w, h, N, "dat/cheby1_approx.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

67
test/src/cheby1_ap_even_odd.c
Wyświetl plik

@ -1,67 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define ORD 5
#define N 1000
#define K 1024
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1];
double Rp = 2.0;
double w[N], mag[N], phi[N], tau[N];
double t[K], h[K];
fft_t pfft;
int k;
int res = cheby1_ap(Rp, ORD-1, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
freqs_resp(b, a, ORD-1, w, N, DSPL_FLAG_LOG|DSPL_FLAG_UNWRAP, mag, phi, tau);
writetxt(w, mag, N, "dat/cheby1_ap_ord4_mag.txt");
writetxt(w, phi, N, "dat/cheby1_ap_ord4_phi.txt");
writetxt(w, tau, N, "dat/cheby1_ap_ord4_tau.txt");
memset(&pfft, 0, sizeof(fft_t));
freqs2time(b, a, ORD-1, 50.0, K, &pfft, t,h);
writetxt(t, h, K, "dat/cheby1_ap_ord4_time.txt");
res = cheby1_ap(Rp, ORD, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD+1; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
freqs_resp(b, a, ORD, w, N, DSPL_FLAG_LOG|DSPL_FLAG_UNWRAP, mag, phi, tau);
writetxt(w, mag, N, "dat/cheby1_ap_ord5_mag.txt");
writetxt(w, phi, N, "dat/cheby1_ap_ord5_phi.txt");
writetxt(w, tau, N, "dat/cheby1_ap_ord5_tau.txt");
memset(&pfft, 0, sizeof(fft_t));
freqs2time(b, a, ORD, 50.0, K, &pfft, t,h);
writetxt(t, h, K, "dat/cheby1_ap_ord5_time.txt");
dspl_free(handle); // free dspl handle
return 0;
}

46
test/src/cheby1_ap_example.c
Wyświetl plik

@ -1,46 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define ORD 4
#define N 1000
#define K 1024
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1];
double Rp = 1.5;
double w[N], mag[N], phi[N], tau[N];
double t[K], h[K];
fft_t pfft;
int k;
int res = cheby1_ap(Rp, ORD, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD+1; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
freqs_resp(b, a, ORD, w, N, DSPL_FLAG_LOG|DSPL_FLAG_UNWRAP, mag, phi, tau);
writetxt(w, mag, N, "dat/cheby1_ap_test_mag.txt");
writetxt(w, phi, N, "dat/cheby1_ap_test_phi.txt");
writetxt(w, tau, N, "dat/cheby1_ap_test_tau.txt");
memset(&pfft, 0, sizeof(fft_t));
freqs2time(b, a, ORD, 50.0, K, &pfft, t,h);
writetxt(t, h, K, "dat/cheby1_ap_test_time.txt");
dspl_free(handle); // free dspl handle
return 0;
}

39
test/src/cheby2_ap_approx.c
Wyświetl plik

@ -1,39 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 1000
#define ORD 5
int main(int argc, char* argv[])
{
double w[N], h[N];
double a[ORD+1], b[ORD+1];
void* handle;
int k;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(0, 2.5, N, DSPL_PERIODIC, w);
cheby2_ap(20.0, ORD, b,a);
freqs_resp(b, a, ORD, w, N, 0, h, NULL, NULL);
for( k =0; k< N; k++)
{
h[k] *= 6.0;
w[k] *= 4.0;
}
writetxt(w, h, N, "dat/cheby2_approx.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

70
test/src/cheby2_ap_even_odd.c
Wyświetl plik

@ -1,70 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define ORD 9
#define N 1000
#define K 1024
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1];
double Rs = 40.0;
double w[N], mag[N], phi[N], tau[N];
double t[K], h[K];
fft_t pfft;
int k;
int res = cheby2_ap(Rs, ORD-1, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
freqs_resp(b, a, ORD-1, w, N,
DSPL_FLAG_LOG|DSPL_FLAG_UNWRAP, mag, phi, tau);
writetxt(w, mag, N, "dat/cheby2_ap_ord6_mag.txt");
writetxt(w, phi, N, "dat/cheby2_ap_ord6_phi.txt");
writetxt(w, tau, N, "dat/cheby2_ap_ord6_tau.txt");
memset(&pfft, 0, sizeof(fft_t));
freqs2time(b, a, ORD-1, 50.0, K, &pfft, t,h);
writetxt(t, h, K, "dat/cheby2_ap_ord6_time.txt");
res = cheby2_ap(Rs, ORD, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD+1; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
freqs_resp(b, a, ORD, w, N,
DSPL_FLAG_LOG|DSPL_FLAG_UNWRAP, mag, phi, tau);
writetxt(w, mag, N, "dat/cheby2_ap_ord7_mag.txt");
writetxt(w, phi, N, "dat/cheby2_ap_ord7_phi.txt");
writetxt(w, tau, N, "dat/cheby2_ap_ord7_tau.txt");
memset(&pfft, 0, sizeof(fft_t));
freqs2time(b, a, ORD, 50.0, K, &pfft, t,h);
writetxt(t, h, K, "dat/cheby2_ap_ord7_time.txt");
dspl_free(handle); // free dspl handle
return 0;
}

46
test/src/cheby2_ap_example.c
Wyświetl plik

@ -1,46 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define ORD 5
#define N 1000
#define K 1024
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1];
double Rs = 50;
double w[N], mag[N], phi[N], tau[N];
double t[K], h[K];
fft_t pfft;
int k;
int res = cheby2_ap(Rs, ORD, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
for(k = 0; k < ORD+1; k++)
printf("b[%2d] = %9.5f a[%2d] = %9.5f\n", k, b[k], k, a[k]);
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
freqs_resp(b, a, ORD, w, N, DSPL_FLAG_LOG|DSPL_FLAG_UNWRAP, mag, phi, tau);
writetxt(w, mag, N, "dat/cheby2_ap_test_mag.txt");
writetxt(w, phi, N, "dat/cheby2_ap_test_phi.txt");
writetxt(w, tau, N, "dat/cheby2_ap_test_tau.txt");
memset(&pfft, 0, sizeof(fft_t));
freqs2time(b, a, ORD, 10.0, K, &pfft, t,h);
writetxt(t, h, K, "dat/cheby2_ap_test_time.txt");
dspl_free(handle); // free dspl handle
return 0;
}

101
test/src/cheby2_ap_zp.c
Wyświetl plik

@ -1,101 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include "dspl.h"
#define ORD 9
#define N 1024
#define SCALE 1.0
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
complex_t z[ORD], p[ORD];
int nz, np, k;
double Rs = 40.0;
printf("\nORD = 9\n");
int res = cheby2_ap_zp(ORD, Rs, z, &nz, p, &np);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
printf("\nChebyshev type 2 zeros:\n");
for(k = 0; k < nz; k++)
printf("(%9.3f, %9.3f)\n", SCALE*RE(z[k]), SCALE*IM(z[k]));
printf("\nChebyshev type 2 poles:\n");
for(k = 0; k < np; k++)
printf("(%9.3f, %9.3f)\n", SCALE*RE(p[k]), SCALE*IM(p[k]));
double alpha[N], sigma[N], omega[N];
double rs, es, beta, shb, chb, ca, sa, den;
char fname[64];
int m;
linspace(0, M_2PI, N, DSPL_SYMMETRIC, alpha);
for(m = 1; m < 12; m++)
{
rs = (double)m*10.0;
es = sqrt(pow(10.0, rs*0.1) - 1.0);
beta = asinh(es)/(double)ORD;
chb = cosh(beta);
shb = sinh(beta);
for(k = 0; k < N; k++)
{
ca = cos(alpha[k]);
sa = sin(alpha[k]);
den = (ca*ca*chb*chb + sa*sa*shb*shb);
sigma[k] = -SCALE * sa * shb/den;
omega[k] = SCALE * ca * chb/den;
}
sprintf(fname, "dat/cheby2_ap_poles_ord9_rs%d.txt", m*10);
writetxt(sigma, omega, N, fname);
}
printf("\nORD = 8\n");
res = cheby2_ap_zp(ORD-1, Rs, z, &nz, p, &np);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
printf("\nChebyshev type 2 zeros:\n");
for(k = 0; k < nz; k++)
printf("(%9.3f, %9.3f)\n", SCALE*RE(z[k]), SCALE*IM(z[k]));
printf("\nChebyshev type 2 poles:\n");
for(k = 0; k < np; k++)
printf("(%9.3f, %9.3f)\n", SCALE*RE(p[k]), SCALE*IM(p[k]));
for(m = 1; m < 12; m++)
{
rs = (double)m*10.0;
es = sqrt(pow(10.0, rs*0.1) - 1.0);
beta = asinh(es)/(double)(ORD-1);
chb = cosh(beta);
shb = sinh(beta);
for(k = 0; k < N; k++)
{
ca = cos(alpha[k]);
sa = sin(alpha[k]);
den = (ca*ca*chb*chb + sa*sa*shb*shb);
sigma[k] = -SCALE * sa * shb/den;
omega[k] = SCALE * ca * chb/den;
}
sprintf(fname, "dat/cheby2_ap_poles_ord8_rs%d.txt", m*10);
writetxt(sigma, omega, N, fname);
}
dspl_free(handle); // free dspl handle
return 0;
}

50
test/src/dft_indexation.c
Wyświetl plik

@ -1,50 +0,0 @@
#include <stdio.h>
#include <string.h>
#include "dspl.h"
#define N 30
#define FS 120
#define F0 -20
int main(int argc, char* argv[])
{
void* handle;
handle = dspl_load();
complex_t s[N]; // Входной сигнал
complex_t X[N]; // ДПФ
double f[N]; // Индексы спектральных отсчетов
double S[N]; // Амплитудный спектр без перестановки
double Ssh[N]; // Амплитудный спектр после перестановки
int n;
// входной сигнал
for(n = 0; n < N; n++)
{
RE(s[n]) = cos(M_2PI * (double)n * F0 / FS);
IM(s[n]) = sin(M_2PI * (double)n * F0 / FS);
}
// ДПФ
dft_cmplx(s, N, X);
// Амплитудный спектр
for(n = 0; n < N; n++)
S[n] = ABS(X[n]);
// Перестановка спектральных отсчетов
fft_shift(S, N, Ssh);
// заполнение массива индексов спектральных отсчетов
linspace(0, N, N, DSPL_PERIODIC, f);
//сохранить данные для построения графика
writetxt(f, S, N, "dat/dft_freq_fig1.txt");
writetxt(f, Ssh, N, "dat/dft_freq_fig5.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

45
test/src/filter_approx.c
Wyświetl plik

@ -1,45 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 1000
int main(int argc, char* argv[])
{
double w[N], r[N], h[N];
double ep2, Gp2;
int ord, k;
void* handle;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(0, 2.5, N, DSPL_PERIODIC, w);
ord = 4;
Gp2 = 0.9;
ep2 = 1.0 / Gp2 -1;
for(k = 0; k < N; k++)
{
r[k] = pow(w[k], (double)(ord));
r[k] *= r[k];
h[k] = 6.0/(1.0 + ep2*r[k]);
w[k] *= 4;
}
writetxt(w, h, N, "dat/butter_approx.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

35
test/src/fourier_series_dirichlet_ex.c
Wyświetl plik

@ -1,35 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 10000
int main(int argc, char* argv[])
{
double t[N], s[N];
void* handle;
int k;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-M_2PI, M_2PI, N, DSPL_PERIODIC, t);
for(k = 0; k < N; k++)
s[k] = sin(1.0/sin(t[k]));
writetxt(t, s, N, "dat/fourier_series_signal0.txt");
linspace(-0.2, 0.2, N, DSPL_PERIODIC, t);
for(k = 0; k < N; k++)
s[k] = sin(1.0/sin(t[k]));
writetxt(t, s, N, "dat/fourier_series_signal1.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

95
test/src/fourier_series_pimp_env.c
Wyświetl plik

@ -1,95 +0,0 @@
#include <stdio.h>
#include <string.h>
#include "dspl.h"
#define N 800
#define K 19
void scale_vector(double* x, int n, double a)
{
int k;
for(k = 0; k < n; k++)
x[k] *= a;
}
void magphi(double* x, int n, double* mag, double* phi)
{
int k;
for(k = 0; k < n/2; k++)
{
mag[k] = fabs(x[k]);
phi[k] = x[k] < 0.0 ? 1.0 : 0.0;
}
for(k = n/2; k < n; k++)
{
mag[k] = fabs(x[k]);
phi[k] = x[k] < 0.0 ? -1.0 : 0.0;
}
}
int main(int argc, char* argv[])
{
double w[N];
double S[N];
double mag[N];
double phi[N];
double wn[K];
double Sn[K];
int n;
void* handle;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-8.0, 8.0, N, DSPL_SYMMETRIC, w);
linspace(-8.0, 8.0, K, DSPL_SYMMETRIC, wn);
sinc(w, N, M_PI*0.5, S);
scale_vector(S, N, 2);
sinc(wn, K, M_PI*0.5, Sn);
scale_vector(Sn, K, 2);
writetxt(w+10, S+10, N - 20, "dat/pimp_env.txt");
writetxt(wn+1, Sn+1, K - 2, "dat/pimp_spec.txt");
magphi(S, N, mag, phi);
writetxt(w+10, mag+10, N - 20, "dat/pimp_mag_env.txt");
writetxt(w+10, mag+10, N - 20, "dat/pimp_mag_env_shift.txt");
writetxt(w+10, phi+10, N - 20, "dat/pimp_phi_env.txt");
magphi(Sn, K, mag, phi);
writetxt(wn+1, mag+1, K-2, "dat/pimp_mag_spec.txt");
writetxt(wn+1, mag+1, K-2, "dat/pimp_mag_spec_shift.txt");
writetxt(wn+1, phi+1, K-2, "dat/pimp_phi_spec.txt");
magphi(S, N, mag, phi);
for(n = 0; n < N; n++)
phi[n] -= w[n]/M_PI;
writetxt(w+10, phi+10, N - 20, "dat/pimp_phi_env_shift.txt");
magphi(Sn, K, mag, phi);
for(n = 0; n < K; n++)
phi[n] -= wn[n]/M_PI;
writetxt(wn+1, phi+1, K - 2, "dat/pimp_phi_spec_shift.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

93
test/src/fourier_series_pimp_q.c
Wyświetl plik

@ -1,93 +0,0 @@
#include <stdio.h>
#include <string.h>
//#include "common.h"
#include "dspl.h"
//#include "plot.h"
#define N 1000
#define T 4.0
#define A 2.0
#define M 41
int main(int argc, char* argv[])
{
double t[N]; // время (сек)
double s[N]; // входной сигнал
complex_t S[M]; // комплексный спектр периодического сигнала
double Smag[M]; // амплитудный спектр периодического сигнала
double w[M]; // частота (рад/c) дискретного спектра
double wc[N]; // частота (рад/с) огибающей спектра
double Sc[N]; // огибающая спектра
double tau; // длительность импульса
// скважность
double Q[3] = {5.0, 2.0, 1.25};
int q, m, n;
char fname[64];
void* handle;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
for(q = 0; q < 3; q++)
{
tau = T/Q[q];
// заполнение массива временных отсчетов
// на 4-x периодах повторения сигнала
// для отображения на осциллограмме
linspace(-T*2.0, T*2.0, N, DSPL_PERIODIC, t);
// 4 периода повторения п-импульса скважности Q[q]
signal_pimp(t, N, A, tau, 0.0, T, s);
// сохранение в текстовый файл
sprintf(fname, "dat/pimp_time_%.2lf.csv", Q[q]);
writetxt(t, s, N, fname);
// заполнение массива временных отсчетов
// на одном периоде повторения сигнала
linspace(-T/2.0, T/2.0, N, DSPL_PERIODIC, t);
// один период повторения п-импульса скважности Q[q]
signal_pimp(t, N, A, tau, 0.0, T, s);
// разложение в ряд Фурье
fourier_series_dec(t, s, N, T, M, w, S);
// Рассчет амплитудного спектра
for(m = 0; m < M; m++)
{
printf("S[%d] = %f %f\n", m, RE(S[m]), IM(S[m]));
Smag[m] = ABS(S[m]);
}
// Сохранение в файл амплитудного спетра для скважности Q[q]
sprintf(fname, "dat/pimp_freq_discrete_%.2lf.csv", Q[q]);
writetxt(w, Smag, M, fname);
// Вектор частот непрерывной огибаюхей вида sin(w/2*tau) / (w/2*T)
linspace(w[0], w[M-1], N, DSPL_SYMMETRIC, wc);
// Расчет огибающей
for(n = 0; n < N; n++)
Sc[n] = (wc[n] == 0.0) ? A/Q[q] : fabs( A * sin(0.5*wc[n]*tau) / (0.5*wc[n] * T));
// сохранение огибающей в файл для скважности Q[q]
sprintf(fname, "dat/pimp_freq_cont_%.2lf.csv", Q[q]);
writetxt(wc, Sc, N, fname);
}
// remember to free the resource
dspl_free(handle);
return 0;
}

58
test/src/fourier_series_pimp_spectrum.c
Wyświetl plik

@ -1,58 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 1000
#define T 4
#define TAU 2
#define K 41
int main(int argc, char* argv[])
{
double t[N]; // время
double s[N]; // исходный сигнал
double w[K]; // массив частоты
complex_t S[K]; // спектр
double mag[K]; // амплитудный спектр
double phi[K]; // фазовый спектр
void* handle;
// test
int n;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
// заполняем массив времени для одного периода сигнала
linspace(-T*0.5, T*0.5, N, DSPL_SYMMETRIC, t);
// сигнал
signal_pimp(t, N, 2.0, TAU, 0.0, T, s);
// расчет комплексного спектра
fourier_series_dec(t, s, N, T, K, w, S);
// Амплитудный и фазовый спектры
for(n = 0; n < K; n++)
{
mag[n] = ABS(S[n]);
phi[n] = atan2(IM(S[n]), RE(S[n]));
}
// Сохранение амплитудного спектра в файл
writetxt(w, mag, K, "dat/fourier_series_pimp_mag.txt");
// Сохранение фазового спектра в файл
writetxt(w, phi, K, "dat/fourier_series_pimp_phi.txt");
// remember to free the resource
dspl_free(handle);
return system("gnuplot -p gnuplot/fourier_series_pimp.plt");
}

118
test/src/fourier_series_prop_mod.c
Wyświetl plik

@ -1,118 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 200000
#define T 4.0
#define A 2.0
#define M 121
#define TAU 1.0
#define DEC 50
int main(int argc, char* argv[])
{
double* t = NULL; // время (сек)
double* s = NULL; // входной сигнал
complex_t* sc = NULL;
complex_t S[M]; // комплексный спектр периодического сигнала
double Smag[M]; // амплитудный спектр периодического сигнала
double Sphi[M]; // фазовый спектр периодического сигнала
double w[M]; // частота (рад/c) дискретного спектра
void* handle;
int k, i;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
t = (double*)malloc(N*sizeof(double));
s = (double*)malloc(N*sizeof(double));
sc = (complex_t*)malloc(N*sizeof(complex_t));
// массив времени
linspace(-T*2.0, T*2.0, N, DSPL_PERIODIC, t);
// 4 периода повторения п-импульса
signal_pimp(t, N, A, TAU, 0.0, T, s);
// расчет спектра
fourier_series_dec(t, s, N, T, M, w, S);
for(k = 0; k < M; k++)
{
Smag[k] = ABS(S[k])/4.0;
Sphi[k] = atan2(IM(S[k]), RE(S[k]));
}
writetxt(w, Smag, M, "dat/fourier_series_prop_spectrum_amp_a.txt");
writetxt(w, Sphi, M, "dat/fourier_series_prop_spectrum_phi_a.txt");
double w0 = 10.0*M_PI;
for(k = 0; k < N; k++)
{
RE(sc[k]) = s[k] * cos(w0*t[k]);
IM(sc[k]) = s[k] * sin(w0*t[k]);
}
// разложение в ряд Фурье
fourier_series_dec_cmplx(t, sc, N, T, M, w, S);
for(k = 0; k < M; k++)
{
Smag[k] = ABS(S[k])/4.0;
Sphi[k] = atan2(IM(S[k]), RE(S[k]));
}
writetxt(w, Smag, M, "dat/fourier_series_prop_spectrum_amp_se.txt");
writetxt(w, Sphi, M, "dat/fourier_series_prop_spectrum_phi_se.txt");
// сохраняю в файл dat/fourier_series_prop_time_a.txt
decimate(t, N, DEC, t, &i);
decimate(s, N, DEC, s, &i);
writetxt(t, s, i, "dat/fourier_series_prop_time_a.txt");
decimate_cmplx(sc, N, DEC, sc, &i);
writetxt_cmplx_re(t, sc, i, "dat/fourier_series_prop_time_se_re.txt");
writetxt_cmplx_im(t, sc, i, "dat/fourier_series_prop_time_se_im.txt");
// массив времени
linspace(-T*2.0, T*2.0, N, DSPL_PERIODIC, t);
// 4 периода повторения п-импульса
signal_pimp(t, N, A, TAU, 0.0, T, s);
for(k = 0; k < N; k++)
{
s[k] *= cos(w0*t[k]);
}
// разложение в ряд Фурье
fourier_series_dec(t, s, N, T, M, w, S);
for(k = 0; k < M; k++)
{
Smag[k] = ABS(S[k])/4.0;
Sphi[k] = atan2(IM(S[k]), RE(S[k]));
}
writetxt(w, Smag, M, "dat/fourier_series_prop_spectrum_amp_sc.txt");
writetxt(w, Sphi, M, "dat/fourier_series_prop_spectrum_phi_sc.txt");
// сохраняю в файл dat/fourier_series_prop_time_sc.txt
decimate(t, N, DEC, t, &i);
decimate(s, N, DEC, s, &i);
writetxt(t, s, i, "dat/fourier_series_prop_time_sc.txt");
free(s);
free(t);
free(sc);
return 0;
}

124
test/src/fourier_series_rec.c
Wyświetl plik

@ -1,124 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
// размер массивов исходных и восстановленных сигналов
#define N 1000
// Период
#define T 2
// Длительность импульса
#define TAU 1
// Максимальное количество коэффициентов ряда Фурье
#define MAX_M 61
// Количество периодов
#define P 3.0
int main(int argc, char* argv[])
{
double t[N]; // время
double s[N]; // исходный сигнал
double x[N]; // восстановленный сигнал
double w[MAX_M];// массив частоты
complex_t S[MAX_M]; // Спектр
complex_t xc[N]; // восстановленный по спектру
// комплексный сигнал
// количество спектральных гармоник усеченного ряда
// Заметим, что 5 гармоник усеченного ряда содержат
// две пары комплексно-сопряженных спектральных компонент
// и одну постоянную составляющую.
// Это означает, что ряд из 5 гармоник в комплексной форме соответствует
// двум значениям a_n и b_n ряда в тригонометрической форме.
// Аналогично M = 21 соответствует 10 значениям a_n и b_n ряда
// в тригонометрической форме.
int M[4] = {5, 9, 21, MAX_M};
void* handle;
int k, n;
char fname[64]; // имя файла
// Загрузка libdspl
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
// Заполняю вектор времени
linspace(-T*0.5*P, T*0.5*P, N, DSPL_PERIODIC, t);
// последовательность прямоугольных импульсов
signal_pimp(t, N, 1.0, TAU, 0.0, T, s);
// сохраняем в файл dat/fourier_series_pimp0.txt
writetxt(t, s, N, "dat/fourier_series_pimp0.txt");
// цикл по разному количеству гармоник усеченного ряда
for(k = 0; k < 4; k++)
{
// расчет спектра для текущего усеченного ряда
fourier_series_dec(t, s, N, T, M[k], w, S);
// нормировка потому что P периодов в исходном сигнале
for(n = 0; n < M[k]; n++)
{
RE(S[n]) /= P;
IM(S[n]) /= P;
}
// восстанавливаю сигнал по усеченному ряду
fourier_series_rec(w, S, M[k], t, N, xc);
// Комплексный восстановленный сигнал имеет очень маленькую
// мнимую часть, поэтому просто берем в вектор x реальную часть
cmplx2re(xc, N, x, NULL);
// сохраняю в файл для последующего построения графика
sprintf(fname, "dat/fourier_series_pimp_rec_%d.txt", M[k]);
writetxt(t, x, N, fname);
}
// Пилообразный сигнал
signal_saw(t, N, 0.5, 0.0, T, s);
for(n = 0; n < N; n++)
s[n] += 0.5;
// сохраняем в файл dat/fourier_series_saw0.txt
writetxt(t, s, N, "dat/fourier_series_saw0.txt");
// цикл по разному количеству гармоник усеченного ряда
for(k = 0; k < 4; k++)
{
// расчет спектра для текущего усеченного ряда
fourier_series_dec(t, s, N, T, M[k], w, S);
// нормировка потому что P периодов в исходном сигнале
for(n = 0; n < M[k]; n++)
{
RE(S[n]) /= P;
IM(S[n]) /= P;
}
// восстанавливаю сигнал по усеченному ряду
fourier_series_rec(w, S, M[k], t, N, xc);
// Комплексный восстановленный сигнал имеет очень маленькую
// мнимую часть, поэтому просто берем в вектор x реальную часть
cmplx2re(xc, N, x, NULL);
// сохраняю в файл для последующего построения графика
sprintf(fname, "dat/fourier_series_saw_rec_%d.txt", M[k]);
writetxt(t, x, N, fname);
}
// remember to free the resource
dspl_free(handle);
return system("gnuplot -p gnuplot/fourier_series_rec.plt");
}

96
test/src/fourier_spectrum_bw.c
Wyświetl plik

@ -1,96 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 1024
#define K 8
#define DT 0.01
#define M 20
int main(int argc, char* argv[])
{
double t[N*K], s[N*K], z[N*K];
double w[N*K], S[N*K], Z[N*K];
double amp[M] = {0.098, -0.160, -0.134, 0.259, 0.248,
-0.525, 0.124, -0.080, 0.179, 0.160,
0.071, -0.055, -0.119, 0.061, 0.297,
0.023, -0.347, 0.085, 0.046, 0.022};
double mu[M] = {478.0, 552.0, 672.0, 189.0, 268.0,
504.0, 824.0, 482.0, 471.0, 643.0,
685.0, 163.0, 870.0, 424.0, 662.0,
611.0, 377.0, 549.0, 824.0, 731.0};
double sigma[M] = {527.0, 1399.0, 1499.0, 543.0, 523.0,
1465.0, 1175.0, 1535.0, 604.0, 1074.0,
890.0, 616.0, 563.0, 1201.0, 1884.0,
1745.0, 1465.0, 1885.0, 1767.0, 1635.0};
complex_t SC[N*K];
int k, m;
double ms;
void* handle;
fft_t pfft;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(0, (double)(N*K), N*K, DSPL_PERIODIC, t);
memset(z, 0, N*K*sizeof(double));
ms = 0.0;
for(k = 0; k < N; k++)
{
for(m =0; m < M; m++)
{
s[k] += amp[m] * exp(-(t[k]-mu[m])*(t[k]-mu[m]) / sigma[m]);
}
s[k] += 2.0 * exp(-(t[k]-(double)N * 0.5)*(t[k]-(double)N * 0.5) / 150.0);
}
for(k = 503; k < 523; k++)
z[k] = 1.8;
for(k = 0; k < N; k++)
t[k] *= DT;
writetxt(t, s, N, "dat/dat_s.txt");
writetxt(t, z, N, "dat/dat_z.txt");
memset(&pfft, 0, sizeof(fft_t));
fft_create(&pfft, N*K);
fft(s, N*K, &pfft, SC);
for(k = 0; k < N*K; k++)
{
S[k] = ABS(SC[k])*DT;
w[k] = -M_PI/DT + M_2PI/DT/(double)(N*K) * (double)k;
}
fft_shift(S, N*K, S);
fft(z, N*K, &pfft, SC);
for(k = 0; k < N*K; k++)
{
Z[k] = ABS(SC[k])*DT;
}
fft_shift(Z, N*K, Z);
writetxt(w, S, N*K, "dat/spec_s.txt");
writetxt(w, Z, N*K, "dat/spec_z.txt");
dspl_free(handle);
return 0;
}

72
test/src/fourier_transform_energy.c
Wyświetl plik

@ -1,72 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 2000
#define A 1
#define TAU 1
int main(int argc, char* argv[])
{
double t[N], w[N], Sr[N], St[N], Se[N], Sg[N], sigma;
void* handle;
int k;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-10, 10, N, DSPL_PERIODIC, t);
linspace(-20*M_PI, 20*M_PI, N, DSPL_PERIODIC, w);
for(k = 0; k < N; k++)
{
Sr[k] = sin(w[k]*TAU*0.5+1E-9) / (w[k]*TAU*0.5+1E-9);
Sr[k] *= Sr[k];
St[k] = sin(w[k]*TAU*0.5+1E-9)/(w[k]*TAU*0.5+1E-9);
St[k] *= St[k];
St[k] *= St[k];
sigma = 0.5;
Sg[k] = w[k]*sigma;
Sg[k] *= w[k]*sigma / 4.0;
Sg[k] = exp(-Sg[k]);
Sg[k] *= Sg[k];
sigma = 2;
Se[k] = sigma*sigma/(w[k]*w[k] + sigma*sigma);
Se[k] *= Se[k];
}
writetxt(w, Sr, N, "dat/fourier_transform_energy_r_lin.txt");
writetxt(w, St, N, "dat/fourier_transform_energy_t_lin.txt");
writetxt(w, Se, N, "dat/fourier_transform_energy_e_lin.txt");
writetxt(w, Sg, N, "dat/fourier_transform_energy_g_lin.txt");
for(k = 0; k < N; k++)
{
Sr[k] = 10.0*log10(Sr[k]);
Sg[k] = 10.0*log10(Sg[k]);
St[k] = 10.0*log10(St[k]);
Se[k] = 10.0*log10(Se[k]);
}
writetxt(w, Sr, N, "dat/fourier_transform_energy_r_log.txt");
writetxt(w, St, N, "dat/fourier_transform_energy_t_log.txt");
writetxt(w, Se, N, "dat/fourier_transform_energy_e_log.txt");
writetxt(w, Sg, N, "dat/fourier_transform_energy_g_log.txt");
// remember to free the resource
dspl_free(handle);
}

114
test/src/fourier_transform_ex.c
Wyświetl plik

@ -1,114 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 2000
#define A 1
int main(int argc, char* argv[])
{
double t[N], w[N], s[N], S[N], PHI[N], sigma;
void* handle;
int k;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-10, 10, N, DSPL_PERIODIC, t);
linspace(-4*M_PI, 4*M_PI, N, DSPL_PERIODIC, w);
sigma = 0.5;
for(k = 0; k < N; k++)
{
s[k] = A * exp(-t[k]*t[k]/(sigma*sigma));
S[k] = A * sigma * sqrt(M_PI) * exp(-w[k]*w[k]*sigma*sigma*0.25);
}
writetxt(t, s, N, "dat/fourier_transform_ex_gauss_time_0.5.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_gauss_freq_0.5.txt");
for(k = 0; k < N; k++)
{
s[k] = A * exp(-fabs(t[k]) * sigma);
S[k] = 2 * A * sigma / (sigma*sigma + w[k]*w[k]);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp2_time_0.5.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp2_freq_0.5.txt");
for(k = 0; k < N; k++)
{
s[k] = t[k] >=0 ? A * exp(-t[k] * sigma) : 0.0;
S[k] = A / sqrt(sigma*sigma + w[k]*w[k]);
PHI[k] = -atan(w[k] / sigma);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp1_time_0.5.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp_freq_mag_0.5.txt");
writetxt(w, PHI, N, "dat/fourier_transform_ex_exp_freq_phi_0.5.txt");
sigma = 1;
for(k = 0; k < N; k++)
{
s[k] = A * exp(-t[k]*t[k]/(sigma*sigma));
S[k] = A * sigma * sqrt(M_PI) * exp(-w[k]*w[k]*sigma*sigma*0.25);
}
writetxt(t, s, N, "dat/fourier_transform_ex_gauss_time_1.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_gauss_freq_1.0.txt");
for(k = 0; k < N; k++)
{
s[k] = A * exp(-fabs(t[k]) * sigma);
S[k] = 2 * A * sigma / (sigma*sigma + w[k]*w[k]);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp2_time_1.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp2_freq_1.0.txt");
for(k = 0; k < N; k++)
{
s[k] = t[k] >=0 ? A * exp(-t[k] * sigma) : 0.0;
S[k] = A / sqrt(sigma*sigma + w[k]*w[k]);
PHI[k] = -atan(w[k] / sigma);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp1_time_1.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp_freq_mag_1.0.txt");
writetxt(w, PHI, N, "dat/fourier_transform_ex_exp_freq_phi_1.0.txt");
sigma = 2;
for(k = 0; k < N; k++)
{
s[k] = A * exp(-t[k]*t[k]/(sigma*sigma));
S[k] = A * sigma * sqrt(M_PI) * exp(-w[k]*w[k]*sigma*sigma*0.25);
}
writetxt(t, s, N, "dat/fourier_transform_ex_gauss_time_2.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_gauss_freq_2.0.txt");
for(k = 0; k < N; k++)
{
s[k] = A * exp(-fabs(t[k]) * sigma);
S[k] = 2 * A * sigma / (sigma*sigma + w[k]*w[k]);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp2_time_2.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp2_freq_2.0.txt");
for(k = 0; k < N; k++)
{
s[k] = t[k] >=0 ? A * exp(-t[k] * sigma) : 0.0;
S[k] = A / sqrt(sigma*sigma + w[k]*w[k]);
PHI[k] = -atan(w[k] / sigma);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp1_time_2.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp_freq_mag_2.0.txt");
writetxt(w, PHI, N, "dat/fourier_transform_ex_exp_freq_phi_2.0.txt");
// remember to free the resource
dspl_free(handle);
return system("gnuplot -p gnuplot/fourier_transform_ex_gauss.plt");;
}

80
test/src/fourier_transform_period_incr.c
Wyświetl plik

@ -1,80 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "dspl.h"
#define N 10000
#define D 10
#define TAU 1
#define K 51
int main(int argc, char* argv[])
{
double *t = NULL, *s = NULL;
double *td = NULL, *sd = NULL;
complex_t S[K];
double w[K], mag[K];
double T;
void* handle;
int n;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
t = (double*)malloc(N*sizeof(double));
s = (double*)malloc(N*sizeof(double));
td = (double*)malloc(N/D*sizeof(double));
sd = (double*)malloc(N/D*sizeof(double));
linspace(-10, 10, N, DSPL_SYMMETRIC, t);
decimate(t, N, D,td, &n);
T = 2.0;
signal_pimp(t, N, 2.0, TAU, 0.0, T, s);
decimate(s, N, D, sd, &n);
writetxt(td, sd, n, "dat/fourier_transform_period_2.0_time.txt");
fourier_series_dec(t, s, N, T, K, w, S);
for(n = 0; n < K; n++)
mag[n] = T*ABS(S[n])/20.0;
writetxt(w, mag, K, "dat/fourier_transform_period_2.0_freq.txt");
T = 4.0;
signal_pimp(t, N, 2.0, TAU, 0.0, T, s);
decimate(s, N, D, sd, &n);
writetxt(td, sd, n, "dat/fourier_transform_period_4.0_time.txt");
fourier_series_dec(t, s, N, T, K, w, S);
for(n = 0; n < K; n++)
mag[n] = T*ABS(S[n])/20.0;
writetxt(w, mag, K, "dat/fourier_transform_period_4.0_freq.txt");
T = 8.0;
signal_pimp(t, N, 2.0, TAU, 0.0, T, s);
decimate(s, N, D, sd, &n);
writetxt(td, sd, n, "dat/fourier_transform_period_8.0_time.txt");
fourier_series_dec(t, s, N, T, K, w, S);
for(n = 0; n < K; n++)
mag[n] = T*ABS(S[n])/20.0;
writetxt(w, mag, K, "dat/fourier_transform_period_8.0_freq.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

55
test/src/fourier_transform_prop_symmetry.c
Wyświetl plik

@ -1,55 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 1000
int main(int argc, char* argv[])
{
double t[N], s[N], S[N], PHI[N], w[N];
complex_t SC[N], sc[N] ;
void* handle;
int k;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-1, 10, N, DSPL_SYMMETRIC, t);
linspace(-2.0*M_2PI, 2.0*M_2PI, N, DSPL_SYMMETRIC, w);
for(k =0; k < N; k++)
{
IM(sc[k]) = 0.0;
if(t[k]<0)
{
s[k] = RE(sc[k]) = 0.0;
}
else
{
s[k] = RE(sc[k]) = exp(-t[k]);
}
}
writetxt(t, s, N/2, "dat/fourier_transform_prop_sym_time.txt");
fourier_integral_cmplx(t, sc, N, N, w, SC);
for(k = 0; k < N; k++)
{
S[k] = ABS(SC[k]);
PHI[k] = ARG(SC[k]);
}
writetxt(w, S, N, "dat/fourier_transform_prop_sym_mag.txt");
writetxt(w, PHI, N, "dat/fourier_transform_prop_sym_phi.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

71
test/src/laplace_tf_3d.c
Wyświetl plik

@ -1,71 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define NW 60
#define NS 30
#define ORD 2
#define R 1.0
#define C 0.5
#define L 2.0
int main(int argc, char* argv[])
{
double w[NW], sigma[NS], habs[NW*NS];
void* handle;
double b[ORD+1] = {1, 0, 0};
double a[ORD+1] = {1, R*C, L*C};
complex_t hs[NW*NS];
complex_t s[NW*NS];
int k, n;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-2, 0, NS, DSPL_SYMMETRIC, sigma);
linspace(-2, 2, NW, DSPL_SYMMETRIC, w);
for(k = 0; k < NW; k++)
{
for(n = 0; n < NS; n++)
{
RE(s[k*NS + n]) = sigma[n];
IM(s[k*NS + n]) = w[k];
}
}
freqs_cmplx(b, a, ORD, s, NW*NS, hs);
//cmplx2re(hs, N, hr, hi);
for(k = 0; k < NW*NS; k++)
{
habs[k] = ABS(hs[k]) > 16.0 ? 16.0 : ABS(hs[k]) ;
}
writetxt_3d(sigma, NS, w, NW, habs, "dat/laplace_tf_3d_abs.txt");
freqs(b,a,ORD, w,NW, hs);
for(k = 0; k<NW; k++)
habs[k] = ABS(hs[k]);
sigma[0] = 0.0;
writetxt_3d(sigma, 1, w, NW, habs, "dat/laplace_tf_3d_hjw.txt");
dspl_free(handle);
return 0;
}

68
test/src/laplace_tf_planes.c
Wyświetl plik

@ -1,68 +0,0 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 200
#define ORD 5
#define SCALEH 2.0
#define SCALEP 0.35
#define SCALES 2.0
int main(int argc, char* argv[])
{
double w[N], sigma[N], hr[N], hi[N];
void* handle;
double b[ORD+1], a[ORD+1];
//complex_t p[ORD], z[ORD];
complex_t hs[N];
int k;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-2.6, 2.6, N, DSPL_SYMMETRIC, w);
butter_ap(2, ORD, b, a);
freqs(b,a, ORD, w, N, hs);
//cmplx2re(hs, N, hr, hi);
for(k = 0; k < N; k++)
{
hr[k] = RE(hs[k]) * SCALES;
hi[k] = IM(hs[k]) * SCALES;
}
memset(sigma, 0, N*sizeof(double));
writetxt_3dline(sigma,w, hr, N, "dat/laplace_tf_planes_3d_re.txt");
writetxt_3dline(sigma,w, hi, N, "dat/laplace_tf_planes_3d_im.txt");
writetxt(w, hr, N, "dat/laplace_tf_planes_2d_re.txt");
writetxt(w, hi, N, "dat/laplace_tf_planes_2d_im.txt");
freqs_resp(b,a, ORD, w, N, DSPL_FLAG_UNWRAP, hr, hi, NULL);
for(k = 0; k < N; k++)
{
hr[k] *= SCALEH;
hi[k] = (hi[k] + M_2PI) * SCALEP;
}
writetxt(w, hr, N, "dat/laplace_tf_planes_2d_abs.txt");
writetxt(w, hi, N, "dat/laplace_tf_planes_2d_phi.txt");
dspl_free(handle);
return 0;
}