kopia lustrzana https://github.com/Dsplib/libdspl-2.0
added many examples and some docs
Dsplib
rodzic
3b27da63d6
commit
b4cb9914ed
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -41,10 +41,13 @@
|
|||
\ref ERROR_CODE_GROUP "код ошибки".<BR>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
\note
|
||||
Для расчета спектра сигнала используетя численное интегрирование
|
||||
Для расчета спектра сигнала используется численное интегрирование
|
||||
исходного сигнала методом трапеций. Данная функция не является
|
||||
вычислительно-эффективной. Для увеличения скорости расчета спектра сигнала
|
||||
эффективной. Для увеличения скорости расчета спектра сигнала
|
||||
целесообразнее использовать алгоритмы дискретного
|
||||
и быстрого преобразования Фурье.
|
||||
<BR>
|
||||
|
|
@ -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$
|
||||
усеченного ряда Фурье.<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] t Указатель на массив временных отсчетов
|
||||
восстановленного сигнала.<BR>
|
||||
Размер вектора `[nt x 1]`.<BR>
|
||||
Память должна быть выделена и заполнена.<BR><BR>
|
||||
<BR><BR>
|
||||
|
||||
\param[in] nt Размер вектора времени и восстановленного сигнала.
|
||||
<BR><BR>
|
||||
|
||||
\param[out] y Указатель на массив восстановленного сигнала.<BR>
|
||||
Размер вектора `[nt x 1]`.<BR>
|
||||
Память должна быть выделена.<BR><BR>
|
||||
|
||||
\return
|
||||
`RES_OK` Массивы нулей и полюсов рассчитаны успешно.<BR>
|
||||
В противном случае
|
||||
\ref ERROR_CODE_GROUP "код ошибки".<BR>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
\note
|
||||
Выходной восстановленный сигнал в общем случае является комплексным.
|
||||
Однако при соблюдении свойств симметрии векторов `w` и `s` относительно
|
||||
нулевой частоты получим мнимую часть элементов вектора `y` на уровне ошибок
|
||||
округления числа с двойной точностью. Ничтожно малую мнимую часть в этом случае
|
||||
можно игнорировать.
|
||||
<BR>
|
||||
|
||||
\author
|
||||
Бахурин Сергей
|
||||
www.dsplib.org
|
||||
|
||||
***************************************************************************** */
|
||||
|
|
|
|||
|
|
@ -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 "код ошибки".<BR>
|
||||
|
||||
|
||||
Пример использования функции `sinc`:
|
||||
|
||||
\include sinc_test.c
|
||||
|
||||
|
||||
<BR><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 произведет построение графика
|
||||
по сохраненным в файлах данным:
|
||||
|
||||
График: `img/sinc_test.png`
|
||||
\image html sinc_test.png
|
||||
|
||||
Скрипт GNUPLOT для построения графиков из текстовых файлов:
|
||||
\include sinc_test.plt
|
||||
|
||||
\author
|
||||
Бахурин Сергей
|
||||
www.dsplib.org
|
||||
|
|
|
|||
|
|
@ -102,3 +102,4 @@ int DSPL_API fourier_series_rec(double* w, complex_t* s, int nw,
|
|||
}
|
||||
return RES_OK;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,21 @@
|
|||
/*!
|
||||
\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 Необходимо выполнить скрипт
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,35 @@
|
|||
|
||||
|
||||
|
||||
|
||||
/*!
|
||||
\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>
|
||||
*/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,24 @@
|
|||
/*! ****************************************************************************
|
||||
\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
|
||||
|
||||
|
||||
|
||||
***************************************************************************** */
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
/*! ****************************************************************************
|
||||
\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
|
||||
|
||||
|
||||
|
||||
***************************************************************************** */
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
/*! ****************************************************************************
|
||||
\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
|
||||
|
||||
***************************************************************************** */
|
||||
|
|
@ -0,0 +1,46 @@
|
|||
#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;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,47 @@
|
|||
#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;
|
||||
}
|
||||
|
|
@ -0,0 +1,67 @@
|
|||
#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;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,46 @@
|
|||
#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;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,39 @@
|
|||
#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;
|
||||
}
|
||||
|
|
@ -0,0 +1,70 @@
|
|||
#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;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,46 @@
|
|||
#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;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,101 @@
|
|||
#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", 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;
|
||||
}
|
||||
|
|
@ -0,0 +1,50 @@
|
|||
#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;
|
||||
}
|
||||
|
|
@ -0,0 +1,128 @@
|
|||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
#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;
|
||||
}
|
||||
|
|
@ -0,0 +1,58 @@
|
|||
#include <stdio.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;
|
||||
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");
|
||||
}
|
||||
|
|
@ -0,0 +1,124 @@
|
|||
#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");
|
||||
}
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
#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;
|
||||
}
|
||||
|
|
@ -0,0 +1,68 @@
|
|||
#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;
|
||||
}
|
||||
Ładowanie…
Reference in New Issue