From b4cb9914ed179b2070892c2e57028a742382f266 Mon Sep 17 00:00:00 2001 From: Dsplib Date: Wed, 3 Oct 2018 23:48:33 +0300 Subject: [PATCH] added many examples and some docs --- dox/doxyfile_ru | 5 +- dox/makedoc.sh | 5 + dspl/dox/ru/fourier_series.dox | 75 ++++++++++- dspl/dox/ru/math.dox | 24 ---- dspl/src/fourier_series.c | 1 + test/bin/gnuplot/fourier_series_pimp.plt | 31 +++++ test/bin/gnuplot/fourier_series_rec.plt | 70 ++++++++++ test/dox/ru/dft_indexation.dox | 21 +++ test/dox/ru/filter_approx.dox | 35 +++++ test/dox/ru/fourier_series_pimp_spectrum.dox | 24 ++++ test/dox/ru/fourier_series_rec.dox | 33 +++++ test/dox/ru/sinc_test.dox | 28 ++++ test/src/butter_ap.c | 46 +++++++ test/src/cheby1_ap_approx.c | 47 +++++++ test/src/cheby1_ap_even_odd.c | 67 ++++++++++ test/src/cheby1_ap_example.c | 46 +++++++ test/src/cheby2_ap_approx.c | 39 ++++++ test/src/cheby2_ap_even_odd.c | 70 ++++++++++ test/src/cheby2_ap_example.c | 46 +++++++ test/src/cheby2_ap_zp.c | 101 +++++++++++++++ test/src/dft_indexation.c | 50 ++++++++ test/src/filter_approx.c | 128 +++++++++++++++++++ test/src/fourier_series_dirichlet_ex.c | 34 +++++ test/src/fourier_series_pimp_spectrum.c | 58 +++++++++ test/src/fourier_series_rec.c | 124 ++++++++++++++++++ test/src/laplace_tf_3d.c | 71 ++++++++++ test/src/laplace_tf_planes.c | 68 ++++++++++ 27 files changed, 1319 insertions(+), 28 deletions(-) create mode 100644 test/bin/gnuplot/fourier_series_pimp.plt create mode 100644 test/bin/gnuplot/fourier_series_rec.plt create mode 100755 test/dox/ru/dft_indexation.dox create mode 100755 test/dox/ru/filter_approx.dox create mode 100644 test/dox/ru/fourier_series_pimp_spectrum.dox create mode 100644 test/dox/ru/fourier_series_rec.dox create mode 100644 test/dox/ru/sinc_test.dox create mode 100755 test/src/butter_ap.c create mode 100755 test/src/cheby1_ap_approx.c create mode 100755 test/src/cheby1_ap_even_odd.c create mode 100755 test/src/cheby1_ap_example.c create mode 100755 test/src/cheby2_ap_approx.c create mode 100755 test/src/cheby2_ap_even_odd.c create mode 100755 test/src/cheby2_ap_example.c create mode 100755 test/src/cheby2_ap_zp.c create mode 100755 test/src/dft_indexation.c create mode 100755 test/src/filter_approx.c create mode 100644 test/src/fourier_series_dirichlet_ex.c create mode 100644 test/src/fourier_series_pimp_spectrum.c create mode 100644 test/src/fourier_series_rec.c create mode 100755 test/src/laplace_tf_3d.c create mode 100755 test/src/laplace_tf_planes.c diff --git a/dox/doxyfile_ru b/dox/doxyfile_ru index 064673b..e349f15 100755 --- a/dox/doxyfile_ru +++ b/dox/doxyfile_ru @@ -435,7 +435,7 @@ LOOKUP_CACHE_SIZE = 0 # normally produced when WARNINGS is set to YES. # The default value is: NO. -EXTRACT_ALL = NO +EXTRACT_ALL = YES # If the EXTRACT_PRIVATE tag is set to YES, all private members of a class will # be included in the documentation. @@ -921,7 +921,8 @@ EXAMPLE_PATH = ../test/src \ # *.h) to filter out the source-files in the directories. If left blank all # files are included. -EXAMPLE_PATTERNS = * +EXAMPLE_PATTERNS = *.c \ + *.dox # If the EXAMPLE_RECURSIVE tag is set to YES then subdirectories will be # searched for input files to be used with the \include or \dontinclude commands diff --git a/dox/makedoc.sh b/dox/makedoc.sh index ce48961..ee0b4a1 100755 --- a/dox/makedoc.sh +++ b/dox/makedoc.sh @@ -1,5 +1,10 @@ #!/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 \; + cd ../ make cd test/bin diff --git a/dspl/dox/ru/fourier_series.dox b/dspl/dox/ru/fourier_series.dox index 9b50145..7b21ba7 100755 --- a/dspl/dox/ru/fourier_series.dox +++ b/dspl/dox/ru/fourier_series.dox @@ -41,10 +41,13 @@ \ref ERROR_CODE_GROUP "код ошибки".
+ + + \note -Для расчета спектра сигнала используетя численное интегрирование +Для расчета спектра сигнала используется численное интегрирование исходного сигнала методом трапеций. Данная функция не является -вычислительно-эффективной. Для увеличения скорости расчета спектра сигнала +эффективной. Для увеличения скорости расчета спектра сигнала целесообразнее использовать алгоритмы дискретного и быстрого преобразования Фурье.
@@ -54,3 +57,71 @@ www.dsplib.org ***************************************************************************** */ + + + + + + +/*! **************************************************************************** +\ingroup DFT_GROUP +\fn int fourier_series_rec(double* w, complex_t* s, int nw, + double *t, int nt, complex_t* y) + +\brief Восстановление сигнала при усечении ряда Фурье + +Функция рассчитывает восстановленный сигнал при усечении ряда Фурье: + +\f[ +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$ + усеченного ряда Фурье.
+ Размер вектора `[nw x 1]`.
+ Память должна быть выделена и заполнена.

+ +\param[in] s Указатель на массив значений спектра + \f$S(\omega_n)\f$.
+ Размер вектора `[nw x 1]`.
+ Память должна быть выделена и заполнена.

+ + +\param[in] nw Количество членов усеченного ряда Фурье.
+ Значение должно быть положительным.

+ +\param[in] t Указатель на массив временных отсчетов + восстановленного сигнала.
+ Размер вектора `[nt x 1]`.
+ Память должна быть выделена и заполнена.

+

+ +\param[in] nt Размер вектора времени и восстановленного сигнала. +

+ +\param[out] y Указатель на массив восстановленного сигнала.
+ Размер вектора `[nt x 1]`.
+ Память должна быть выделена.

+ +\return + `RES_OK` Массивы нулей и полюсов рассчитаны успешно.
+ В противном случае + \ref ERROR_CODE_GROUP "код ошибки".
+ + + + + +\note +Выходной восстановленный сигнал в общем случае является комплексным. +Однако при соблюдении свойств симметрии векторов `w` и `s` относительно +нулевой частоты получим мнимую часть элементов вектора `y` на уровне ошибок +округления числа с двойной точностью. Ничтожно малую мнимую часть в этом случае +можно игнорировать. +
+ +\author + Бахурин Сергей + www.dsplib.org + +***************************************************************************** */ diff --git a/dspl/dox/ru/math.dox b/dspl/dox/ru/math.dox index 0bd0496..c893aa4 100644 --- a/dspl/dox/ru/math.dox +++ b/dspl/dox/ru/math.dox @@ -1,4 +1,3 @@ - /*! **************************************************************************** \ingroup SPEC_MATH_COMMON_GROUP \fn int sinc(double* x, int n, double a, double* y) @@ -27,29 +26,6 @@ \ref ERROR_CODE_GROUP "код ошибки".
-Пример использования функции `sinc`: - -\include sinc_test.c - - -

-В каталоге `dat` будут созданы три файла:
- -
-sinc_test_1.0.txt	  параметр a = 1.0
-sinc_test_pi.txt	  параметр a = pi
-sinc_test_2pi.txt	  параметр a = 2*pi
-
- -Кроме того программа GNUPLOT произведет построение графика -по сохраненным в файлах данным: - -График: `img/sinc_test.png` -\image html sinc_test.png - -Скрипт GNUPLOT для построения графиков из текстовых файлов: -\include sinc_test.plt - \author Бахурин Сергей www.dsplib.org diff --git a/dspl/src/fourier_series.c b/dspl/src/fourier_series.c index 03bc02e..bb715d6 100755 --- a/dspl/src/fourier_series.c +++ b/dspl/src/fourier_series.c @@ -102,3 +102,4 @@ int DSPL_API fourier_series_rec(double* w, complex_t* s, int nw, } return RES_OK; } + diff --git a/test/bin/gnuplot/fourier_series_pimp.plt b/test/bin/gnuplot/fourier_series_pimp.plt new file mode 100644 index 0000000..b7d13d1 --- /dev/null +++ b/test/bin/gnuplot/fourier_series_pimp.plt @@ -0,0 +1,31 @@ + +set terminal wxt 0 size 560,480 enhanced font 'Verdana,8' position 0,0 +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 + + +set terminal pngcairo size 560,480 enhanced font 'Verdana,8' +set output 'img/fourier_series_pimp.png' +set multiplot layout 2,1 rowsfirst +# +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 +# +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 diff --git a/test/bin/gnuplot/fourier_series_rec.plt b/test/bin/gnuplot/fourier_series_rec.plt new file mode 100644 index 0000000..272d199 --- /dev/null +++ b/test/bin/gnuplot/fourier_series_rec.plt @@ -0,0 +1,70 @@ +### Start multiplot (2x2 layout) + +set terminal wxt 0 size 800,640 enhanced font 'Verdana,8' position 0,0 +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 + + + + +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 diff --git a/test/dox/ru/dft_indexation.dox b/test/dox/ru/dft_indexation.dox new file mode 100755 index 0000000..e1f1b39 --- /dev/null +++ b/test/dox/ru/dft_indexation.dox @@ -0,0 +1,21 @@ +/*! + \example dft_indexation.c + \brief Расчет данных для построения рисунка 1 и рисунка 5 раздела + + "Индексация и перестановка спектральных отсчетов дискретного преобразования Фурье" + . +
+ Программа рассчитает \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`.
+ + Данные расчета сохраняются в текстовые файлы для построения графиков с + использованием пакетов tikz и pgfplots системы LaTeX.
+ + Для построения графиков при помощи пакета GNUPlot Необходимо выполнить скрипт + +*/ diff --git a/test/dox/ru/filter_approx.dox b/test/dox/ru/filter_approx.dox new file mode 100755 index 0000000..b008c74 --- /dev/null +++ b/test/dox/ru/filter_approx.dox @@ -0,0 +1,35 @@ + + + + +/*! + \example filter_approx.c + \brief Расчет данных для построения графиков раздела + « + Расчет аналогового фильтра. Постановка задачи и способы аппроксимации АЧХ нормированного ФНЧ + ». + + + Программа сохраняет текстовые файлы данных для построения графиков аппроксимирующих функций \f$ F_N^2(\omega)\f$ и + квадрата АЧХ \f$ |H(j\omega)|^2\f$ нормированных аналоговых фильтров нижних частот: + +
+
`dat/butter_r.txt`
+
+ Файл значений квадрата аппроксимирующей функции фильтра Баттерворта +
+
`dat/butter_h.txt`
+
+ Файл значений квадрата АЧХ фильтра Баттерворта в линейном масштабе +
+
`dat/butter_hdb.txt`
+
+ Файл значений квадрата АЧХ фильтра Баттерворта в логарифмическом масштабе (дБ) +
+
+*/ + + + + + diff --git a/test/dox/ru/fourier_series_pimp_spectrum.dox b/test/dox/ru/fourier_series_pimp_spectrum.dox new file mode 100644 index 0000000..8988019 --- /dev/null +++ b/test/dox/ru/fourier_series_pimp_spectrum.dox @@ -0,0 +1,24 @@ +/*! **************************************************************************** +\example fourier_series_pimp_spectrum.c +Расчет амплитудного и фазового спектра периодической последовательности +прямоугольных импульсов +
+ +Программа формирует один период последовательности прямоугольных импульсов +и производит расчет амплитудного и фазового спектра. + +Результаты расчета сохраняются в файлы и выводятся на график программой GNUPlot + + +Скрипт GNUPLOT `fourier_series_pimp.plt` +для построения графиков из текстовых файлов: + +\include fourier_series_pimp.plt + + +График будет отображен на экране и сохранен в файл `img/fourier_series_pimp.png` +\image html fourier_series_pimp.png + + + +***************************************************************************** */ diff --git a/test/dox/ru/fourier_series_rec.dox b/test/dox/ru/fourier_series_rec.dox new file mode 100644 index 0000000..9227267 --- /dev/null +++ b/test/dox/ru/fourier_series_rec.dox @@ -0,0 +1,33 @@ +/*! **************************************************************************** +\example fourier_series_rec.c +Программа производит расчет данных для представления +периодической последовательности прямоугольных импульсов и пилообразного сигнала +усеченным рядом Фурье. +
+ +Программа формирует `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 + + + +***************************************************************************** */ diff --git a/test/dox/ru/sinc_test.dox b/test/dox/ru/sinc_test.dox new file mode 100644 index 0000000..3473901 --- /dev/null +++ b/test/dox/ru/sinc_test.dox @@ -0,0 +1,28 @@ +/*! **************************************************************************** +\example sinc_test.c Пример использования функции \ref sinc + +Программа производит расчет функции +\f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$ +для 3 различных значений параметра \f$ a\f$. +
+ +В каталоге `dat` будут созданы три файла:
+ +
+sinc_test_1.0.txt	  параметр a = 1.0
+sinc_test_pi.txt	  параметр a = pi
+sinc_test_2pi.txt	  параметр a = 2*pi
+
+ + +Кроме того программа GNUPLOT произведет построение графика +по сохраненным в файлах данным. + +Скрипт GNUPLOT `sinc_test.plt` для построения графиков из текстовых файлов: +\include sinc_test.plt + + +График будет отображен на экране и сохранен в файл `img/sinc_test.png` +\image html sinc_test.png + +***************************************************************************** */ diff --git a/test/src/butter_ap.c b/test/src/butter_ap.c new file mode 100755 index 0000000..9142757 --- /dev/null +++ b/test/src/butter_ap.c @@ -0,0 +1,46 @@ +#include +#include +#include +#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; +} + + diff --git a/test/src/cheby1_ap_approx.c b/test/src/cheby1_ap_approx.c new file mode 100755 index 0000000..7a3fc9d --- /dev/null +++ b/test/src/cheby1_ap_approx.c @@ -0,0 +1,47 @@ +#include +#include +#include +#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; +} diff --git a/test/src/cheby1_ap_even_odd.c b/test/src/cheby1_ap_even_odd.c new file mode 100755 index 0000000..37a548c --- /dev/null +++ b/test/src/cheby1_ap_even_odd.c @@ -0,0 +1,67 @@ +#include +#include +#include +#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; +} + + diff --git a/test/src/cheby1_ap_example.c b/test/src/cheby1_ap_example.c new file mode 100755 index 0000000..4c118bf --- /dev/null +++ b/test/src/cheby1_ap_example.c @@ -0,0 +1,46 @@ +#include +#include +#include +#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; +} + + diff --git a/test/src/cheby2_ap_approx.c b/test/src/cheby2_ap_approx.c new file mode 100755 index 0000000..ba0159b --- /dev/null +++ b/test/src/cheby2_ap_approx.c @@ -0,0 +1,39 @@ +#include +#include +#include +#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; +} diff --git a/test/src/cheby2_ap_even_odd.c b/test/src/cheby2_ap_even_odd.c new file mode 100755 index 0000000..67d499f --- /dev/null +++ b/test/src/cheby2_ap_even_odd.c @@ -0,0 +1,70 @@ +#include +#include +#include +#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; +} + + diff --git a/test/src/cheby2_ap_example.c b/test/src/cheby2_ap_example.c new file mode 100755 index 0000000..f35d9c3 --- /dev/null +++ b/test/src/cheby2_ap_example.c @@ -0,0 +1,46 @@ +#include +#include +#include +#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; +} + + diff --git a/test/src/cheby2_ap_zp.c b/test/src/cheby2_ap_zp.c new file mode 100755 index 0000000..0b4793a --- /dev/null +++ b/test/src/cheby2_ap_zp.c @@ -0,0 +1,101 @@ +#include +#include +#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", k, 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", k, 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", k, 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", k, 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; +} diff --git a/test/src/dft_indexation.c b/test/src/dft_indexation.c new file mode 100755 index 0000000..d606946 --- /dev/null +++ b/test/src/dft_indexation.c @@ -0,0 +1,50 @@ +#include +#include +#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; +} diff --git a/test/src/filter_approx.c b/test/src/filter_approx.c new file mode 100755 index 0000000..d84224b --- /dev/null +++ b/test/src/filter_approx.c @@ -0,0 +1,128 @@ +#include +#include +#include +#include "dspl.h" + +#define N 8000 + + + +void r2hdb(double* r, double* h, double* hdb, double ep2, int n) +{ + int k; + for(k = 0; k < n; k++) + { + r[k] *= r[k]; + h[k] = 1.0/(1.0 + ep2 * r[k]); + hdb[k] = 10.0*log10(h[k]); + } +} + + +int main(int argc, char* argv[]) +{ + double w[N]; // время (сек) + double r[N]; // входной сигнал + double h[N], hdb[N]; + double ep2, Rp, Rs, es2, m; + int ord, k; + + void* handle; // DSPL handle + handle = dspl_load(); // Load DSPL function + + // заполняю массив частот в логарифмическом формате от 0.01 до 100 + logspace(-2, 2, N, DSPL_PERIODIC, w); + + + ord = 4; // порядок фильтра + Rp = 0.3; // Неравномерность в полосе пропускания (дБ) + Rs = 15.0; // Уровень подавления в полосе заграждения (дБ) + + // Параметры ep^2 и es^2 + ep2 = pow(10.0, Rp*0.1) - 1.0; + es2 = pow(10.0, Rs*0.1) - 1.0; + + // вывод на печать параметров ep^2 и es^2 + printf("ep^2 = %.4f\n", ep2); + printf("es^2 = %.4f\n", es2); + + + + + //********************************************************************** + // Расчет F_N^2(w) и |H(jw)|^2 фильтра Баттерворта + //********************************************************************** + for(k = 0; k < N; k++) + { + r[k] = pow(w[k], (double)(ord)); + } + r2hdb(r, h, hdb, ep2, N); + + // сохранение в файлы результатов расчета фильтра Баттерворта + writetxt(w, r, N, "dat/butter_r.txt"); + writetxt(w, h, N, "dat/butter_h.txt"); + writetxt(w, hdb, N, "dat/butter_hdb.txt"); + + + + + + //********************************************************************** + // Расчет F_N^2(w) и |H(jw)|^2 фильтра Чебышева 1-го рода + //********************************************************************** + cheby_poly1(w, N, ord, r); + r2hdb(r, h, hdb, ep2, N); + // сохранение в файлы результатов расчета фильтра Чебышева 1-го рода + writetxt(w, r, N, "dat/cheby1_r.txt"); + writetxt(w, h, N, "dat/cheby1_h.txt"); + writetxt(w, hdb, N, "dat/cheby1_hdb.txt"); + + + + + + //********************************************************************** + // Расчет F_N^2(w) и |H(jw)|^2 фильтра Чебышева 2-го рода + //********************************************************************** + for(k = 0; k < N; k++) + { + w[k] = 1.0 / w[k]; + } + cheby_poly1(w, N, ord, r); + for(k = 0; k < N; k++) + { + r[k] =1.0 / (r[k] *r[k]); + h[k] = 1.0/(1.0 + es2*r[k]); + hdb[k] = 10.0*log10(h[k]); + } + logspace(-2, 2, N, DSPL_PERIODIC, w); + // сохранение в файлы результатов расчета фильтра Чебышева 1-го рода + writetxt(w, r, N, "dat/cheby2_r.txt"); + writetxt(w, h, N, "dat/cheby2_h.txt"); + writetxt(w, hdb, N, "dat/cheby2_hdb.txt"); + + + + + + //********************************************************************** + // Расчет F_N^2(w) и |H(jw)|^2 эллиптического фильтра + //********************************************************************** + ellip_modulareq(Rp, Rs, ord, &m); // пересчет эллиптического модуля + printf("modular m = %.3f\n", m); // вывод на печать + // расчет эллиптической рациональной функции + ellip_rat(w, N, ord, m, r); + r2hdb(r, h, hdb, ep2, N); + // сохранение в файлы результатов расчета эллиптического фильтра + writetxt(w, r, N, "dat/ellip_r.txt"); + writetxt(w, h, N, "dat/ellip_h.txt"); + writetxt(w, hdb, N, "dat/ellip_hdb.txt"); + + + + dspl_free(handle); // free dspl handle + return 0; +} + + + diff --git a/test/src/fourier_series_dirichlet_ex.c b/test/src/fourier_series_dirichlet_ex.c new file mode 100644 index 0000000..f87fc59 --- /dev/null +++ b/test/src/fourier_series_dirichlet_ex.c @@ -0,0 +1,34 @@ +#include +#include +#include +#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; +} diff --git a/test/src/fourier_series_pimp_spectrum.c b/test/src/fourier_series_pimp_spectrum.c new file mode 100644 index 0000000..1bc10c3 --- /dev/null +++ b/test/src/fourier_series_pimp_spectrum.c @@ -0,0 +1,58 @@ +#include +#include +#include +#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; + int k, 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"); +} diff --git a/test/src/fourier_series_rec.c b/test/src/fourier_series_rec.c new file mode 100644 index 0000000..95394e2 --- /dev/null +++ b/test/src/fourier_series_rec.c @@ -0,0 +1,124 @@ +#include +#include +#include +#include +#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"); +} diff --git a/test/src/laplace_tf_3d.c b/test/src/laplace_tf_3d.c new file mode 100755 index 0000000..8ca4716 --- /dev/null +++ b/test/src/laplace_tf_3d.c @@ -0,0 +1,71 @@ +#include +#include +#include +#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 +#include +#include +#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; +}