diff --git a/_release/dspl.c b/_release/dspl.c new file mode 100644 index 0000000..875ecaf --- /dev/null +++ b/_release/dspl.c @@ -0,0 +1,602 @@ +/* +* Copyright (c) 2015-2020 Sergey Bakhurin +* Digital Signal Processing Library [http://dsplib.org] +* +* This file is part of libdspl-2.0. +* +* is free software: you can redistribute it and/or modify +* it under the terms of the GNU Lesser General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* DSPL is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU Lesser General Public License +* along with Foobar. If not, see . +*/ + + + +#ifdef WIN_OS +#include +#endif /* WIN_OS */ + +#ifdef LINUX_OS +#include +#endif /* LINUX_OS */ + + +#include +#include "dspl.h" + + +#ifndef BUILD_LIB + +p_acos_cmplx acos_cmplx ; +p_addlog addlog ; +p_array_scale_lin array_scale_lin ; +p_asin_cmplx asin_cmplx ; + +p_butter_ap butter_ap ; +p_bessel_i0 bessel_i0 ; +p_bilinear bilinear ; +p_butter_ap_zp butter_ap_zp ; + +p_cheby_poly1 cheby_poly1 ; +p_cheby_poly2 cheby_poly2 ; +p_cheby1_ap cheby1_ap ; +p_cheby1_ap_zp cheby1_ap_zp ; +p_cheby2_ap cheby2_ap ; +p_cheby2_ap_wp1 cheby2_ap_wp1 ; +p_cheby2_ap_zp cheby2_ap_zp ; +p_cmplx2re cmplx2re ; +p_concat concat ; +p_conv conv ; +p_conv_cmplx conv_cmplx ; +p_conv_fft conv_fft ; +p_conv_fft_cmplx conv_fft_cmplx ; +p_cos_cmplx cos_cmplx ; + +p_decimate decimate ; +p_decimate_cmplx decimate_cmplx ; +p_dft dft ; +p_dft_cmplx dft_cmplx ; +p_dmod dmod ; +p_dspl_info dspl_info ; + +p_ellip_acd ellip_acd ; +p_ellip_acd_cmplx ellip_acd_cmplx ; +p_ellip_ap ellip_ap ; +p_ellip_ap_zp ellip_ap_zp ; +p_ellip_asn ellip_asn ; +p_ellip_asn_cmplx ellip_asn_cmplx ; +p_ellip_cd ellip_cd ; +p_ellip_cd_cmplx ellip_cd_cmplx ; +p_ellip_landen ellip_landen ; +p_ellip_modulareq ellip_modulareq ; +p_ellip_rat ellip_rat ; +p_ellip_sn ellip_sn ; +p_ellip_sn_cmplx ellip_sn_cmplx ; + +p_farrow_lagrange farrow_lagrange ; +p_farrow_spline farrow_spline ; +p_fft fft ; +p_fft_abs fft_abs ; +p_fft_abs_cmplx fft_abs_cmplx ; +p_fft_cmplx fft_cmplx ; +p_fft_create fft_create ; +p_fft_free fft_free ; +p_fft_mag fft_mag ; +p_fft_mag_cmplx fft_mag_cmplx ; +p_fft_shift fft_shift ; +p_fft_shift_cmplx fft_shift_cmplx ; +p_filter_freq_resp filter_freq_resp ; +p_filter_iir filter_iir ; +p_filter_ws1 filter_ws1 ; +p_filter_zp2ab filter_zp2ab ; +p_find_max_abs find_max_abs ; +p_fir_linphase fir_linphase ; +p_flipip flipip ; +p_flipip_cmplx flipip_cmplx ; +p_fourier_integral_cmplx fourier_integral_cmplx ; +p_fourier_series_dec fourier_series_dec ; +p_fourier_series_dec_cmplx fourier_series_dec_cmplx ; +p_fourier_series_rec fourier_series_rec ; +p_freqs freqs ; +p_freqs_cmplx freqs_cmplx ; +p_freqs2time freqs2time ; +p_freqz freqz ; + +p_gnuplot_close gnuplot_close ; +p_gnuplot_cmd gnuplot_cmd ; +p_gnuplot_create gnuplot_create ; +p_gnuplot_open gnuplot_open ; +p_goertzel goertzel ; +p_goertzel_cmplx goertzel_cmplx ; +p_group_delay group_delay ; + +p_histogram histogram ; +p_histogram_norm histogram_norm ; + +p_idft_cmplx idft_cmplx ; +p_ifft_cmplx ifft_cmplx ; +p_iir iir ; + +p_linspace linspace ; +p_log_cmplx log_cmplx ; +p_logspace logspace ; +p_low2bp low2bp ; +p_low2bs low2bs ; +p_low2high low2high ; +p_low2low low2low ; + +p_matrix_eig_cmplx matrix_eig_cmplx ; +p_matrix_eye matrix_eye ; +p_matrix_eye_cmplx matrix_eye_cmplx ; +p_matrix_mul matrix_mul ; +p_matrix_print matrix_print ; +p_matrix_print_cmplx matrix_print_cmplx ; +p_matrix_transpose matrix_transpose ; +p_matrix_transpose_cmplx matrix_transpose_cmplx ; +p_matrix_transpose_hermite matrix_transpose_hermite ; +p_mean mean ; +p_mean_cmplx mean_cmplx ; +p_minmax minmax ; + +p_ones ones ; + +p_phase_delay phase_delay ; +p_poly_z2a_cmplx poly_z2a_cmplx ; +p_polyroots polyroots ; +p_polyval polyval ; +p_polyval_cmplx polyval_cmplx ; +p_psd_bartlett psd_bartlett ; +p_psd_bartlett_cmplx psd_bartlett_cmplx ; +p_psd_periodogram psd_periodogram ; +p_psd_periodogram_cmplx psd_periodogram_cmplx ; +p_psd_welch psd_welch ; +p_psd_welch_cmplx psd_welch_cmplx ; + +p_randb randb ; +p_randb2 randb2 ; +p_randi randi ; +p_randn randn ; +p_randn_cmplx randn_cmplx ; +p_random_init random_init ; +p_randu randu ; +p_ratcompos ratcompos ; +p_re2cmplx re2cmplx ; +p_readbin readbin ; + +p_signal_pimp signal_pimp ; +p_signal_saw signal_saw ; +p_sin_cmplx sin_cmplx ; +p_sinc sinc ; +p_sine_int sine_int ; +p_sqrt_cmplx sqrt_cmplx ; +p_std std ; +p_std_cmplx std_cmplx ; + +p_trapint trapint ; +p_trapint_cmplx trapint_cmplx ; + +p_unwrap unwrap ; + +p_vector_dot vector_dot ; +p_verif verif ; +p_verif_data_gen verif_data_gen ; +p_verif_cmplx verif_cmplx ; +p_verif_str verif_str ; +p_verif_str_cmplx verif_str_cmplx ; + +p_window window ; +p_writebin writebin ; +p_writetxt writetxt ; +p_writetxt_3d writetxt_3d ; +p_writetxt_3dline writetxt_3dline ; +p_writetxt_cmplx writetxt_cmplx ; +p_writetxt_cmplx_im writetxt_cmplx_im ; +p_writetxt_cmplx_re writetxt_cmplx_re ; +p_writetxt_int writetxt_int ; + +p_xcorr xcorr ; +p_xcorr_cmplx xcorr_cmplx ; + + +#ifdef WIN_OS +#define LOAD_FUNC(fn) \ + fname = #fn;\ + fn = (p_##fn)GetProcAddress(handle, fname);\ + if(! fn) goto exit_label; +#endif + + + +#ifdef LINUX_OS +#define LOAD_FUNC(fn) \ + fname = #fn;\ + fn = (p_##fn)dlsym(handle, fname);\ + if ((error = dlerror()) != NULL) goto exit_label +#endif + + + +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup SYS_LOADING_GROUP +\fn void* dspl_load() +\brief Perform dynamic linking and load libdspl-2.0 functions. + +This function attempts to link to the library `libdspl.dll` in +Windows system and the `libdspl.so` library on the Linux system. +The library is assumed to be in the same directory as the application. +user, or the path to the library is registered in the operating path variables +system. + +Upon successful binding and loading of library functions, the handle is returned +libraries, as well as in the address space of the application appear +pointers to libdspl-2.0 functions. + +\note +The returned handle is of type `void *`, which can be cast on Windows +to type `HINSTANCE`. In practice, this is not necessary, because this +the type is cast to `HINSTANCE` automatically if the compiler flag is set, +indicating that the application is being built on Windows. + +An example of a simple program that implements dynamic binding with DSPL-2.0. + +\code +#include +#include +#include "dspl.h" + +int main (int argc, char* argv[]) +{ + void * hdspl; // DSPL handle + hdspl = dspl_load (); // Dynamic linking + + // Check the pointer. If `NULL`, then the link failed + if (! hdspl) + { + printf ("libdspl loading error! \n"); + return -1; + } + + // The link was successful, you can call the functions of DSPL-2.0 + + //Before correctly terminating the application, you must unlink + //library and clear memory. + dspl_free(hdspl); + + return 0; +} +\endcode + +\author Bakhurin Sergey. www.dsplib.org +***************************************************************************** */ +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup SYS_LOADING_GROUP +\fn void* dspl_load() +\brief Произвести динамическую линковку и загрузить функции libdspl-2.0. + +Данная функция производит попытку связывания с библиотекой `libdspl.dll` в +системе Windows и с библиотекой `libdspl.so` в системе Linux. +Предполагается, что библиотека находится в одной директории с приложением +пользователя, или путь к библиотеке прописан в переменных пути операционной +системы. + +При удачном связывании и загрузке функций библиотеки возвращается хэндл +библиотеки, а также в адресном пространстве приложения появляются +указатели на функции libdspl-2.0. + +\note +Возвращаемый хэндл имеет тип `void*`, который в ОС Windows может быть приведен +к типу `HINSTANCE`. На практике необходимости в этом, нет, потому что данный +тип приводится к `HINSTANCE` автоматически, если выставлен флаг компилятора, +указывающий, что сборка приложения производится в ОС Windows. + +Пример простейшей программы реализующей динамическое связывание с DSPL-2.0. + +\code +#include +#include +#include "dspl.h" + +int main(int argc, char* argv[]) +{ + void* hdspl; // DSPL хэндл + hdspl = dspl_load(); // Динамическая линковка + + // Проверяем указатель. Если `NULLL`, то линковка прошла неудачно + if(!hdspl) + { + printf("libdspl loading error!\n"); + return -1; + } + + // Линковка прошла успешно можно вызывать функции DSPL-2.0 + + // Перед корректным завершением приложения необходимо разлинковать + // библиотеку и очистить память. + dspl_free(hdspl); + + return 0; +} +\endcode + +\author Бахурин Сергей. www.dsplib.org +***************************************************************************** */ +#endif +void* dspl_load() +{ + char* fname; + #ifdef WIN_OS + HINSTANCE handle; + handle = LoadLibrary(TEXT("libdspl.dll")); + if (!handle) + { + printf("libdspl.dll loading ERROR!\n"); + return NULL; + } + #endif /* WIN_OS */ + + + #ifdef LINUX_OS + char* error; + void *handle; + /* open the *.so */ + handle = dlopen ("./libdspl.so", RTLD_LAZY); + if (!handle) + { + printf("libdspl.so loading ERROR!\n"); + return NULL; + } + #endif /* LINUX_OS */ + + LOAD_FUNC(acos_cmplx); + LOAD_FUNC(addlog); + LOAD_FUNC(array_scale_lin); + LOAD_FUNC(asin_cmplx); + + LOAD_FUNC(bessel_i0); + LOAD_FUNC(bilinear); + LOAD_FUNC(butter_ap); + LOAD_FUNC(butter_ap_zp); + + LOAD_FUNC(cheby_poly1); + LOAD_FUNC(cheby_poly2); + LOAD_FUNC(cheby1_ap); + LOAD_FUNC(cheby1_ap_zp); + LOAD_FUNC(cheby2_ap); + LOAD_FUNC(cheby2_ap_wp1); + LOAD_FUNC(cheby2_ap_zp); + LOAD_FUNC(cmplx2re); + LOAD_FUNC(concat); + LOAD_FUNC(conv); + LOAD_FUNC(conv_cmplx); + LOAD_FUNC(conv_fft); + LOAD_FUNC(conv_fft_cmplx); + LOAD_FUNC(cos_cmplx); + + LOAD_FUNC(decimate); + LOAD_FUNC(decimate_cmplx); + LOAD_FUNC(dft); + LOAD_FUNC(dft_cmplx); + LOAD_FUNC(dmod); + LOAD_FUNC(dspl_info); + + LOAD_FUNC(ellip_acd); + LOAD_FUNC(ellip_acd_cmplx); + LOAD_FUNC(ellip_ap); + LOAD_FUNC(ellip_ap_zp); + LOAD_FUNC(ellip_asn); + LOAD_FUNC(ellip_asn_cmplx); + LOAD_FUNC(ellip_cd); + LOAD_FUNC(ellip_cd_cmplx); + LOAD_FUNC(ellip_landen); + LOAD_FUNC(ellip_modulareq); + LOAD_FUNC(ellip_rat); + LOAD_FUNC(ellip_sn); + LOAD_FUNC(ellip_sn_cmplx); + + LOAD_FUNC(farrow_lagrange); + LOAD_FUNC(farrow_spline); + LOAD_FUNC(fft); + LOAD_FUNC(fft_cmplx); + LOAD_FUNC(fft_create); + LOAD_FUNC(fft_free); + LOAD_FUNC(fft_mag); + LOAD_FUNC(fft_mag_cmplx); + LOAD_FUNC(fft_shift); + LOAD_FUNC(fft_shift_cmplx); + LOAD_FUNC(filter_freq_resp); + LOAD_FUNC(filter_iir); + LOAD_FUNC(filter_ws1); + LOAD_FUNC(filter_zp2ab); + LOAD_FUNC(find_max_abs); + LOAD_FUNC(fir_linphase); + LOAD_FUNC(flipip); + LOAD_FUNC(flipip_cmplx); + LOAD_FUNC(fourier_integral_cmplx); + LOAD_FUNC(fourier_series_dec); + LOAD_FUNC(fourier_series_dec_cmplx); + LOAD_FUNC(fourier_series_rec); + LOAD_FUNC(freqz); + LOAD_FUNC(freqs); + LOAD_FUNC(freqs_cmplx); + LOAD_FUNC(freqs2time); + + LOAD_FUNC(gnuplot_close); + LOAD_FUNC(gnuplot_cmd); + LOAD_FUNC(gnuplot_create); + LOAD_FUNC(gnuplot_open); + LOAD_FUNC(goertzel); + LOAD_FUNC(goertzel_cmplx); + LOAD_FUNC(group_delay); + + LOAD_FUNC(histogram); + LOAD_FUNC(histogram_norm); + + LOAD_FUNC(idft_cmplx); + LOAD_FUNC(ifft_cmplx); + LOAD_FUNC(iir); + + LOAD_FUNC(linspace); + LOAD_FUNC(log_cmplx); + LOAD_FUNC(logspace); + LOAD_FUNC(low2bp); + LOAD_FUNC(low2bs); + LOAD_FUNC(low2high); + LOAD_FUNC(low2low); + + LOAD_FUNC(matrix_eig_cmplx); + LOAD_FUNC(matrix_eye); + LOAD_FUNC(matrix_eye_cmplx); + LOAD_FUNC(matrix_mul); + LOAD_FUNC(matrix_print); + LOAD_FUNC(matrix_print_cmplx); + LOAD_FUNC(matrix_transpose); + LOAD_FUNC(matrix_transpose_cmplx); + LOAD_FUNC(matrix_transpose_hermite); + LOAD_FUNC(mean); + LOAD_FUNC(mean_cmplx); + LOAD_FUNC(minmax); + + LOAD_FUNC(ones); + + LOAD_FUNC(phase_delay); + LOAD_FUNC(poly_z2a_cmplx); + LOAD_FUNC(polyroots); + LOAD_FUNC(polyval); + LOAD_FUNC(polyval_cmplx); + LOAD_FUNC(psd_bartlett); + LOAD_FUNC(psd_bartlett_cmplx); + LOAD_FUNC(psd_periodogram); + LOAD_FUNC(psd_periodogram_cmplx); + LOAD_FUNC(psd_welch); + LOAD_FUNC(psd_welch_cmplx); + + LOAD_FUNC(randi); + LOAD_FUNC(randb); + LOAD_FUNC(randb2); + LOAD_FUNC(randn); + LOAD_FUNC(randn_cmplx); + LOAD_FUNC(random_init); + LOAD_FUNC(randu); + LOAD_FUNC(ratcompos); + LOAD_FUNC(re2cmplx); + LOAD_FUNC(readbin); + + LOAD_FUNC(signal_pimp); + LOAD_FUNC(signal_saw); + LOAD_FUNC(sin_cmplx); + LOAD_FUNC(sinc); + LOAD_FUNC(sine_int); + LOAD_FUNC(sqrt_cmplx); + LOAD_FUNC(std); + LOAD_FUNC(std_cmplx); + + LOAD_FUNC(trapint); + LOAD_FUNC(trapint_cmplx); + + LOAD_FUNC(unwrap); + + LOAD_FUNC(vector_dot); + LOAD_FUNC(verif); + LOAD_FUNC(verif_data_gen); + LOAD_FUNC(verif_cmplx); + LOAD_FUNC(verif_str); + LOAD_FUNC(verif_str_cmplx); + + LOAD_FUNC(window); + LOAD_FUNC(writebin); + LOAD_FUNC(writetxt); + LOAD_FUNC(writetxt_3d); + LOAD_FUNC(writetxt_3dline); + LOAD_FUNC(writetxt_cmplx); + LOAD_FUNC(writetxt_cmplx_im); + LOAD_FUNC(writetxt_cmplx_re); + LOAD_FUNC(writetxt_int); + + LOAD_FUNC(xcorr); + LOAD_FUNC(xcorr_cmplx); + + + #ifdef WIN_OS + return (void*)handle; + exit_label: + printf("function %s loading ERROR!\n", fname); + if(handle) + FreeLibrary(handle); + return NULL; + #endif /* WIN_OS */ + + + #ifdef LINUX_OS + return handle; + exit_label: + printf("function %s loading ERROR!\n", fname); + if(handle) + dlclose(handle); + return NULL; + #endif /* LINUX_OS */ +} + + + + + + +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup SYS_LOADING_GROUP +\fn void dspl_free(void* handle) +\brief Cleans up the previously linked DSPL-2.0 dynamic library. + +This cross-platform function clears the library `libdspl.dll` in +Windows system and from the library `libdspl.so` on the Linux system. +After cleaning the library, all functions will become unavailable. + +\param [in] handle +Handle of the previously linked DSPL-2.0 library. \n +This pointer can be `NULL`, in this case no action +are being produced. + +\author Bakhurin Sergey. www.dsplib.org +***************************************************************************** */ +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup SYS_LOADING_GROUP +\fn void dspl_free(void* handle) +\brief Очищает связанную ранее динамическую библиотеку DSPL-2.0. + +Данная кроссплатформенная функция производит очистку библиотеки `libdspl.dll` в +системе Windows и с библиотеки `libdspl.so` в системе Linux. +После очистки библиотеки все функции станут недоступны. + +\param[in] handle +Хэндл прилинкованной ранее библиотеки DSPL-2.0. \n +Данный указатель может быть `NULL`, в этом случае никакие действия не +производятся. + +\author Бахурин Сергей. www.dsplib.org +**************************************************************************** */ +#endif +void dspl_free(void* handle) +{ + #ifdef WIN_OS + FreeLibrary((HINSTANCE)handle); + #endif /* WIN_OS */ + + #ifdef LINUX_OS + dlclose(handle); + #endif /* LINUX_OS */ +} + +#endif /* BUILD_LIB */ diff --git a/_release/dspl.h b/_release/dspl.h new file mode 100644 index 0000000..b5b7bef --- /dev/null +++ b/_release/dspl.h @@ -0,0 +1,1572 @@ +/* +* Copyright (c) 2015-2020 Sergey Bakhurin +* Digital Signal Processing Library [http://dsplib.org] +* +* This file is part of libdspl-2.0. +* +* is free software: you can redistribute it and/or modify +* it under the terms of the GNU Lesser General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* DSPL is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU Lesser General Public License +* along with Foobar. If not, see . +*/ + + +#ifndef DSPL_H +#define DSPL_H + +#include + +/* math const definition */ +#ifndef M_PI + #define M_PI 3.1415926535897932384626433832795 +#endif + +#ifndef M_2PI + #define M_2PI 6.283185307179586476925286766559 +#endif + + +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup TYPES_GROUP +\typedef complex_t +\brief Complex data type. + +libdspl-2.0 describes complex numbers data type as an array +of two `double` elements. +First element sets real part, second --- imaginary part. + +For example: + +\code{.cpp} + complex_t z; + z[0] = 1.0; + z[1] = -2.0; +\endcode + +Variable `z = 1-2j`, here `j` - imaginary unit. + +For the convenience of working with complex numbers implemented +special macros: \ref RE, \ref IM, \ref ABSSQR +***************************************************************************** */ +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup TYPES_GROUP +\typedef complex_t +\brief . + + libdspl-2.0 + `double`. + + , - . + +: + +\code{.cpp} + complex_t z; + z[0] = 1.0; + z[1] = -2.0; +\endcode + + `z = 1-2j`, `j` - . + + + : \ref RE, \ref IM, \ref ABSSQR +***************************************************************************** */ +#endif +typedef double complex_t[2]; + + + +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup DFT_GROUP +\struct fft_t +\brief Fast Fourier Transform Object Data Structure + +The structure stores pointers to twiddle factors and arrays of intermediate +data of the fast Fourier transform algorithm. + +The libdspl-2.0 library uses an FFT algorithm for composite size. + +\param n +The size of the FFT vector for which memory is allocated +in the structure arrays. \n +The parameter `n` must be equal to an integer power of two (radix 2). \n \n + +\param w +Pointer to the vector of twiddle factors. \n +The size of the vector is `[n x 1]`. \n +The memory must be allocated and an array of twiddle factors +must be filled with the \ref fft_create function. \n\n + +\param t0 +Pointer to the vector of intermediate results of the FFT algorithm. \n +The size of the vector is `[n x 1]`. \n +Memory must be allocated by \ref fft_create function. \n\n + +\param t1 +Pointer to the vector of intermediate results. \n +The size of the vector is `[n x 1]`. \n +The memory must be allocated with the \ref fft_create function. \n\n +The structure is populated with the \ref fft_create function once +before using the FFT algorithm. \n +A pointer to an object of this structure may be +reused when calling FFT functions. \n +Before exiting the program, dedicated memory for twiddle factors and arrays of +intermediate data must be cleared by the \ref fft_free function. + +For example: + +\code +fft_t pfft = {0}; // Structure fft_t and clear all fields +int n = 64; // FFT size + +int err; + +// Create and fill FFT structure for 64-points FFT +err = fft_create(&pfft, n); + +// FFT calculation here +// FFT calculation here one more +// ... + +// Clear fft structure +fft_free(&pfft); +\endcode + +\note +It is important to note that if the object `fft_t` was created for the FFT size +equal to` n`, it can only be used for FFT of size `n`. \n \n +Its also worth noting that the FFT functions independently control the size, +and independently allocate the memory of the FFT object, if necessary. +So if you call any function using the `fft_t` structure with filled +data for the FFT length `k` for calculating the FFT of length`n`, +then the structure arrays will be automatically recreated for the length `n`. + +\author Sergey Bakhurin www.dsplib.org +***************************************************************************** */ +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup DFT_GROUP +\struct fft_t +\brief + + + . + + libdspl-2.0 + +\param n + , . \n + `n` . \n \n + +\param w + . \n + `[n x 1]`. \n + + \ref fft_create. \n \n + +\param t0 + . \n + `[n x 1]`. \n + \ref fft_create. \n \n + +\param t1 + . \n + `[n x 1]`. \n + \ref fft_create. \n \n + \ref fft_create + . \n + + . \n + + + \ref fft_free. : +\code +fft_t pfft = {0}; // fft_t +int n = 64; // +int err; + +// 64- +err = fft_create(&pfft, n); + +// +// +// ... + +// +fft_free(&pfft); +\endcode + +\note + , `fft_t` `n`, + `n`. \n\n + , , + . + `fft_t` + `k` `n`, + `n`. + +\author + . +www.dsplib.org +***************************************************************************** */ +#endif +typedef struct +{ + complex_t* w; + complex_t* t0; + complex_t* t1; + int n; +} fft_t; + + + +#define RAND_TYPE_MRG32K3A 0x00000001 +#define RAND_TYPE_MT19937 0x00000002 +#define RAND_MT19937_NN 312 + +#ifdef DOXYGEN_ENGLISH + +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup SPEC_MATH_RAND_GEN_GROUP +\struct random_t +\brief . + + + . : +\li MRG32K3A -- 32 [1]. +\li MT19937-64 -- 64- + + + [2, 3]. + +\note +[1] Pierre L'Ecuyer, (1999) Good Parameters and Implementations for Combined + Multiple Recursive Random Number Generators. Operations Research + 47(1):159-164. https://doi.org/10.1287/opre.47.1.159 \n\n +[2] T. Nishimura, ``Tables of 64-bit Mersenne Twisters // ACM Transactions + on Modeling and Computer Simulation 10. (2000) 348--357. \n\n +[3] M. Matsumoto and T. Nishimura Mersenne Twister: a 623-dimensionally + equidistributed uniform pseudorandom number generator // ACM Transactions + on Modeling and Computer Simulation 8. (Jan. 1998) 3--30. \n\n + +\param mrg32k3a_seed + MRG32K3A. \n \n + +\param mrg32k3a_x + MRG32K3A. \n \n + +\param mrg32k3a_y + MRG32K3A. \n \n + +\param mt19937_mt + MT19937-64. \n \n + +\param mt19937_mti + MT19937-64. \n \n + + `random_init` + . + +\author . www.dsplib.org +***************************************************************************** */ +#endif +typedef struct +{ + + double mrg32k3a_seed; + double mrg32k3a_x[3]; + double mrg32k3a_y[3]; + + /* The array for the MT19937 state vector */ + unsigned long long mt19937_mt[RAND_MT19937_NN]; + int mt19937_mti; + + int type; + +}random_t; + + + +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup TYPES_GROUP +\def RE(x) +\brief Macro sets real part of the complex number. + +Example: +\code{.cpp} + complex_t z; + RE(z) = 1.0; + IM(z) = -2.0; +\endcode + +Variable `z = 1-2j`, here `j` - imaginary unit. + +This macro can be used to return +real part of the complex number: + +\code{.cpp} + complex_t z = {3.0, -4.0}; + double r; + r = RE(z); +\endcode +In this example `z = 3-4i`, +but variable `r` will keep 3. +***************************************************************************** */ +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup TYPES_GROUP +\def RE(x) +\brief . + +: +\code{.cpp} + complex_t z; + RE(z) = 1.0; + IM(z) = -2.0; +\endcode + + `z = 1-2j`, `j` --- . + +, + : + +\code{.cpp} + complex_t z = {3.0, -4.0}; + double r; + r = RE(z); +\endcode + `z = 3-4i`, `r` + 3. +***************************************************************************** */ +#endif +#define RE(x) (x[0]) + + + +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup TYPES_GROUP +\def IM(x) +\brief Macro sets imaginary part of the complex number. + +Example: +\code{.cpp} + complex_t z; + RE(z) = 1.0; + IM(z) = -2.0; +\endcode + +Variable `z = 1-2j`, here `j` - imaginary unit. + +This macro can be used to return +imaginary part of the complex number: +\code{.cpp} + complex_t z = {3.0, -4.0}; + double r; + r = IM(z); +\endcode +In this example `z = 3-4i`, +but variable `r` will keep -4. +***************************************************************************** */ + + +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup TYPES_GROUP +\def IM(x) +\brief . + +: +\code{.cpp} + complex_t z; + RE(z) = 1.0; + IM(z) = -2.0; +\endcode + + `z = 1-2j`, `j` - . + +, + : +\code{.cpp} + complex_t z = {3.0, -4.0}; + double r; + r = IM(z); +\endcode + `z = 3-4i`, + `r` -4. +***************************************************************************** */ +#endif +#define IM(x) (x[1]) + + + +#define SQR(x) ((x) * (x)) + + +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup TYPES_GROUP +\def ABSSQR(x) +\brief +The macro returns the square of the modulus of a complex number `x`. + +Square of the modulus of a complex number \f$ x = a + j b \f$ equals: + +\f[ + |x|^2 = x x^* = a^2 + b^2. +\f] + +Example: +\code{.cpp} + complex_t z; + double y; + RE(z) = 1.0; + IM(z) = -2.0; + y = ABSSQR(z); +\endcode + +Variable `z = 1-2j`, here `j` - imaginary unit, but variable `y = 5`. +***************************************************************************** */ +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup TYPES_GROUP +\def ABSSQR(x) +\brief `x`. + + \f$ x = a + j b \f$ : + +\f[ + |x|^2 = x x^* = a^2 + b^2. +\f] + +: +\code{.cpp} + complex_t z; + double y; + RE(z) = 1.0; + IM(z) = -2.0; + y = ABSSQR(z); +\endcode + + `z = 1-2j`, `j` - , `y = 5`. +***************************************************************************** */ +#endif +#define ABSSQR(x) ((SQR(RE(x))) + (SQR(IM(x)))) + + + + +#define ABS(x) sqrt((ABSSQR(x))) + + +#define ARG(x) atan2(IM(x), RE(x)) + + +#define CMRE(a,b) ((RE(a)) * (RE(b)) - (IM(a)) * (IM(b))) + + +#define CMIM(a,b) ((RE(a)) * (IM(b)) + (IM(a)) * (RE(b))) + + +#define CMCONJRE(a, b) ((RE(a)) * (RE(b)) + (IM(a)) * (IM(b))) + + +#define CMCONJIM(a, b) ((IM(a)) * (RE(b)) - (RE(a)) * (IM(b))) + + + +#define RES_OK 0 + +/* Error codes */ +/* A 0x01xxxxxx*/ +#define ERROR_ARG_PARAM 0x01180716 +/* B 0x02xxxxxx*/ +/* C 0x03xxxxxx*/ +/* D 0x04xxxxxx*/ +#define ERROR_DAT_TYPE 0x04012020 +#define ERROR_DIV_ZERO 0x04102226 +/* E 0x05xxxxxx*/ +#define ERROR_ELLIP_MODULE 0x05121315 +/* F 0x06xxxxxx*/ +#define ERROR_FFT_SIZE 0x06062021 +#define ERROR_FILTER_A0 0x06090100 +#define ERROR_FILTER_APPROX 0x06090116 +#define ERROR_FILTER_FT 0x06090620 +#define ERROR_FILTER_ORD 0x06091518 +#define ERROR_FILTER_ORD_BP 0x06091519 +#define ERROR_FILTER_RP 0x06091816 +#define ERROR_FILTER_RS 0x06091819 +#define ERROR_FILTER_TYPE 0x06092025 +#define ERROR_FILTER_WP 0x06092316 +#define ERROR_FILTER_WS 0x06092319 +#define ERROR_FNAME 0x06140113 +#define ERROR_FOPEN 0x06151605 +#define ERROR_FREAD_SIZE 0x06180501 +#define ERROR_FS 0x06190000 +#define ERROR_FWRITE_SIZE 0x06231820 +/* G 0x07xxxxxx*/ +#define ERROR_GNUPLOT_CREATE 0x07161203 +#define ERROR_GNUPLOT_FNPNG 0x07161206 +#define ERROR_GNUPLOT_TERM 0x07161220 +/* H 0x08xxxxxx*/ +/* I 0x09xxxxxx*/ +#define ERROR_INF 0x09140600 +/* J 0x10xxxxxx*/ +/* K 0x11xxxxxx*/ +/* L 0x12xxxxxx*/ +#define ERROR_LAPACK 0x12011601 +/* M 0x13xxxxxx*/ +#define ERROR_MALLOC 0x13011212 +#define ERROR_MATRIX_SIZE 0x13011926 +#define ERROR_MIN_MAX 0x13091413 +/* N 0x14xxxxxx*/ +#define ERROR_NAN 0x14011400 +#define ERROR_NEGATIVE 0x14050701 +/* O 0x15xxxxxx*/ +#define ERROR_OVERLAP 0x15220412 +/* P 0x16xxxxxx*/ +#define ERROR_POLY_AN 0x16150114 +#define ERROR_POLY_ORD 0x16151518 +#define ERROR_PTR 0x16201800 +/* Q 0x17xxxxxx*/ +/* R 0x18xxxxxx*/ +#define ERROR_RAND_SIGMA 0x18011909 +#define ERROR_RAND_TYPE 0x18012009 +#define ERROR_RESAMPLE_RATIO 0x18051801 +#define ERROR_RESAMPLE_FRAC_DELAY 0x18050604 +/* S 0x19xxxxxx*/ +#define ERROR_SIZE 0x19092605 +#define ERROR_SYM_TYPE 0x19251320 +/* T 0x20xxxxxx*/ +/* U 0x21xxxxxx*/ +#define ERROR_UNWRAP 0x21142318 +/* V 0x22xxxxxx*/ +/* W 0x23xxxxxx*/ +#define ERROR_WIN_PARAM 0x23091601 +#define ERROR_WIN_SYM 0x23091925 +#define ERROR_WIN_TYPE 0x23092025 +/* X 0x24xxxxxx*/ +#define ERROR_XCORR_FLAG 0x24031518 +/* Y 0x25xxxxxx*/ +/* Z 0x26xxxxxx*/ + +#define DAT_MASK 0x00000001 +#define DAT_DOUBLE 0x00000000 +#define DAT_COMPLEX 0x00000001 + +#define DSPL_MATRIX_BLOCK 32 + + +#define DSPL_SYMMETRIC 0x00000000 +#define DSPL_PERIODIC 0x00000001 + +#define DSPL_FLAG_DIGITAL 0x00000000 +#define DSPL_FLAG_ANALOG 0x00000001 +#define DSPL_FLAG_LOGMAG 0x00000002 +#define DSPL_FLAG_UNWRAP 0x00000004 +#define DSPL_FLAG_FFT_SHIFT 0x00000008 +#define DSPL_FLAG_PSD_TWOSIDED DSPL_FLAG_FFT_SHIFT + + + + +#define DSPL_WIN_SYM_MASK 0x00000001 +#define DSPL_WIN_MASK 0x00FFFFFE + +#define DSPL_WIN_SYMMETRIC DSPL_SYMMETRIC +#define DSPL_WIN_PERIODIC DSPL_PERIODIC + + +#define DSPL_WIN_BARTLETT 0x00000004 +#define DSPL_WIN_BARTLETT_HANN 0x00000008 +#define DSPL_WIN_BLACKMAN 0x00000010 +#define DSPL_WIN_BLACKMAN_HARRIS 0x00000040 +#define DSPL_WIN_BLACKMAN_NUTTALL 0x00000080 +#define DSPL_WIN_FLAT_TOP 0x00000100 +#define DSPL_WIN_GAUSSIAN 0x00000400 +#define DSPL_WIN_HAMMING 0x00000800 +#define DSPL_WIN_HANN 0x00001000 +#define DSPL_WIN_LANCZOS 0x00004000 +#define DSPL_WIN_NUTTALL 0x00008000 +#define DSPL_WIN_RECT 0x00010000 +#define DSPL_WIN_COS 0x00040000 +#define DSPL_WIN_CHEBY 0x00080000 +#define DSPL_WIN_KAISER 0x00100000 + + +#define DSPL_FILTER_TYPE_MASK 0x000000FF +#define DSPL_FILTER_LPF 0x00000001 +#define DSPL_FILTER_HPF 0x00000002 +#define DSPL_FILTER_BPASS 0x00000004 +#define DSPL_FILTER_BSTOP 0x00000008 + +#define DSPL_FILTER_APPROX_MASK 0x0000FF00 +#define DSPL_FILTER_BUTTER 0x00000100 +#define DSPL_FILTER_CHEBY1 0x00000200 +#define DSPL_FILTER_CHEBY2 0x00000400 +#define DSPL_FILTER_ELLIP 0x00000800 + + +#define DSPL_XCORR_NOSCALE 0x00000000 +#define DSPL_XCORR_BIASED 0x00000001 +#define DSPL_XCORR_UNBIASED 0x00000002 + + + +#define ELLIP_ITER 16 +#define ELLIP_MAX_ORD 24 + +#define DSPL_VERIF_FAILED 1 +#define DSPL_VERIF_SUCCESS 0 + +#define PLOT_HOLD 0x00000001 + +#define VERIF_STR_BUF 128 +#define VERIF_STR_LEN 48 +#define VERIF_CHAR_POINT 46 +#define VERIF_LEVEL_COMPLEX 1E-11 +#define VERIF_LEVEL_DOUBLE 1E-12 + + + +#ifdef __cplusplus + extern "C" { +#endif + + + +#ifdef BUILD_LIB + /* Declare DSPL_API for Windows OS */ + #ifdef WIN_OS + #define DSPL_API __declspec(dllexport) + #endif /* WIN_OS */ + /* Declare DSPL_API for LINUX OS */ + #ifdef LINUX_OS + #define DSPL_API + #endif /* LINUX_OS */ +#endif /* BUILD_DLL */ + +#define COMMA , + + +#ifdef BUILD_LIB + #define DECLARE_FUNC(type, fn, param)\ + type DSPL_API fn(param); +#endif + +#ifndef BUILD_LIB + #define DECLARE_FUNC( type, fn, param)\ + typedef type (*p_##fn)(param);\ + extern p_##fn fn; + +#endif + + + +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, acos_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, addlog, char* str + COMMA char* fn); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, array_scale_lin, double* x + COMMA int n + COMMA double xmin + COMMA double xmax + COMMA double dx + COMMA double h + COMMA double* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, asin_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, bessel_i0, double* x + COMMA int n + COMMA double* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, bilinear, double* bs + COMMA double* as + COMMA int ord + COMMA double* bz + COMMA double* az); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, butter_ap, double + COMMA int + COMMA double* + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, butter_ap_zp, int + COMMA double + COMMA complex_t* + COMMA int* + COMMA complex_t* + COMMA int*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cheby_poly1, double* + COMMA int + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cheby_poly2, double* + COMMA int + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cheby1_ap, double + COMMA int + COMMA double* + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cheby1_ap_zp, int + COMMA double + COMMA complex_t* + COMMA int* + COMMA complex_t* + COMMA int*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cheby2_ap, double rs + COMMA int ord + COMMA double* b + COMMA double* a); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cheby2_ap_wp1, double rp + COMMA double rs + COMMA int ord + COMMA double* b + COMMA double* a); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cheby2_ap_zp, int + COMMA double + COMMA complex_t* + COMMA int* + COMMA complex_t* + COMMA int*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cmplx2re, complex_t* + COMMA int + COMMA double* + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, concat, void* + COMMA size_t + COMMA void* + COMMA size_t + COMMA void*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, conv, double* + COMMA int + COMMA double* + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, conv_cmplx, complex_t* + COMMA int + COMMA complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, conv_fft, double* a + COMMA int na + COMMA double* b + COMMA int nb + COMMA fft_t* pfft + COMMA int nfft + COMMA double* c); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, conv_fft_cmplx, complex_t* a + COMMA int na + COMMA complex_t* b + COMMA int nb + COMMA fft_t* pfft + COMMA int nfft + COMMA complex_t* c); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, cos_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, decimate, double* x + COMMA int n + COMMA int d + COMMA double* y + COMMA int* cnt); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, decimate_cmplx, complex_t* x + COMMA int n + COMMA int d + COMMA complex_t* y + COMMA int* cnt); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, dft, double* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, dft_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(double, dmod, double + COMMA double); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(void, dspl_info, void); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_acd, double* w + COMMA int n + COMMA double k + COMMA double* u); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_acd_cmplx, complex_t* w + COMMA int n + COMMA double k + COMMA complex_t* u); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_ap, double rp + COMMA double rs + COMMA int ord + COMMA double* b + COMMA double* a); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_ap_zp, int ord + COMMA double rp + COMMA double rs + COMMA complex_t* z + COMMA int* nz + COMMA complex_t* p + COMMA int* np); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_asn, double* w + COMMA int n + COMMA double k + COMMA double* u); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_asn_cmplx, complex_t* w + COMMA int n + COMMA double k + COMMA complex_t* u); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_cd, double* u + COMMA int n + COMMA double k + COMMA double* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_cd_cmplx, complex_t* u + COMMA int n + COMMA double k + COMMA complex_t* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_landen, double k + COMMA int n + COMMA double* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_modulareq, double rp + COMMA double rs + COMMA int ord + COMMA double* k); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_rat, double* w + COMMA int n + COMMA int ord + COMMA double k + COMMA double* u); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_sn, double* u + COMMA int n + COMMA double k + COMMA double* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ellip_sn_cmplx, complex_t* u + COMMA int n + COMMA double k + COMMA complex_t* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, farrow_lagrange, double* + COMMA int + COMMA double + COMMA double + COMMA double + COMMA double** + COMMA int*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, farrow_spline, double* + COMMA int + COMMA double + COMMA double + COMMA double + COMMA double** + COMMA int*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft, double* + COMMA int + COMMA fft_t* + COMMA complex_t* ); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_abs, double* x + COMMA int n + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* mag + COMMA double* freq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_abs_cmplx, complex_t* x + COMMA int n + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* mag + COMMA double* freq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_cmplx, complex_t* + COMMA int + COMMA fft_t* + COMMA complex_t* ); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_create, fft_t* + COMMA int); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(void, fft_free, fft_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_mag, double* x + COMMA int n + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* mag + COMMA double* freq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_mag_cmplx, complex_t* x + COMMA int n + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* mag + COMMA double* freq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_shift, double* + COMMA int n + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fft_shift_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, filter_freq_resp, double* b + COMMA double* a + COMMA int ord + COMMA double* w + COMMA int n + COMMA int flag + COMMA double* mag + COMMA double* phi + COMMA double* tau); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, filter_iir, double* + COMMA double* + COMMA int + COMMA double* + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(double, filter_ws1, int ord + COMMA double rp + COMMA double rs + COMMA int type); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, filter_zp2ab, complex_t* + COMMA int + COMMA complex_t* + COMMA int + COMMA int + COMMA double* + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, find_max_abs, double* a + COMMA int n + COMMA double* m + COMMA int* ind); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fir_linphase, int ord + COMMA double w0 + COMMA double w1 + COMMA int filter_type + COMMA int wintype + COMMA double winparam + COMMA double* h); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, flipip, double* + COMMA int); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, flipip_cmplx, complex_t* + COMMA int); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fourier_integral_cmplx, double* t + COMMA complex_t* s + COMMA int nt + COMMA int nw + COMMA double* w + COMMA complex_t* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fourier_series_dec, double* + COMMA double* + COMMA int + COMMA double + COMMA int + COMMA double* + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fourier_series_dec_cmplx, double* t + COMMA complex_t* s + COMMA int nt + COMMA double period + COMMA int nw + COMMA double* w + COMMA complex_t* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, fourier_series_rec, double* + COMMA complex_t* + COMMA int + COMMA double* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, freqs, double* + COMMA double* + COMMA int + COMMA double* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, freqs_cmplx, double* b + COMMA double* a + COMMA int ord + COMMA complex_t* s + COMMA int n + COMMA complex_t* h); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, freqs2time, double* + COMMA double* + COMMA int + COMMA double + COMMA int + COMMA fft_t* + COMMA double* + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, freqz, double* + COMMA double* + COMMA int + COMMA double* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(void, gnuplot_close, void* h); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(void, gnuplot_cmd, void* h + COMMA char* cmd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, gnuplot_create, int argc + COMMA char* argv[] + COMMA int w + COMMA int h + COMMA char* fn_png + COMMA void** hplot); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, gnuplot_open, void** hplot); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, goertzel, double* + COMMA int + COMMA int* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, goertzel_cmplx, complex_t* + COMMA int + COMMA int* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, group_delay, double* b + COMMA double* a + COMMA int ord + COMMA int flag + COMMA double* w + COMMA int n + COMMA double* tau); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, histogram, double* x + COMMA int n + COMMA int nh + COMMA double* pedges + COMMA double* ph); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, histogram_norm, double* y + COMMA int n + COMMA int nh + COMMA double* x + COMMA double* w); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, idft_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ifft_cmplx, complex_t* + COMMA int + COMMA fft_t* + COMMA complex_t* ); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, iir, double rp + COMMA double rs + COMMA int ord + COMMA double w0 + COMMA double w1 + COMMA int type + COMMA double* b + COMMA double* a); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, linspace, double + COMMA double + COMMA int + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, log_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, logspace, double + COMMA double + COMMA int + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, low2bp, double* b + COMMA double* a + COMMA int ord + COMMA double w0 + COMMA double wpl + COMMA double wph + COMMA double* beta + COMMA double* alpha); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, low2bs, double* b + COMMA double* a + COMMA int ord + COMMA double w0 + COMMA double wsl + COMMA double wsh + COMMA double* beta + COMMA double* alpha); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, low2high, double* b + COMMA double* a + COMMA int ord + COMMA double w0 + COMMA double w1 + COMMA double* beta + COMMA double* alpha); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, low2low, double* b + COMMA double* a + COMMA int ord + COMMA double w0 + COMMA double w1 + COMMA double* beta + COMMA double* alpha); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_eig_cmplx, complex_t* a + COMMA int n + COMMA complex_t* v + COMMA int* info); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_eye, double* a + COMMA int n + COMMA int m); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_eye_cmplx, complex_t* a + COMMA int n + COMMA int m); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_mul, double* a + COMMA int na + COMMA int ma + COMMA double* b + COMMA int nb + COMMA int mb + COMMA double* c); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_print, double* a + COMMA int n + COMMA int m + COMMA const char* name + COMMA const char* format); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_print_cmplx, complex_t* a + COMMA int n + COMMA int m + COMMA const char* name + COMMA const char* format); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_transpose, double* a + COMMA int n + COMMA int m + COMMA double* b); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_transpose_cmplx, complex_t* a + COMMA int n + COMMA int m + COMMA complex_t* b); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, matrix_transpose_hermite, complex_t* a + COMMA int n + COMMA int m + COMMA complex_t* b); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, mean, double* x + COMMA int n + COMMA double* m); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, mean_cmplx, complex_t* x + COMMA int n + COMMA complex_t* m); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, minmax, double* x + COMMA int n + COMMA double* xmin + COMMA double* xmax); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ones, double* x + COMMA int n); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, phase_delay, double* b + COMMA double* a + COMMA int ord + COMMA int flag + COMMA double* w + COMMA int n + COMMA double* tau); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, poly_z2a_cmplx, complex_t* + COMMA int + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, polyroots, double* a + COMMA int ord + COMMA complex_t* r + COMMA int* info); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, polyval, double* + COMMA int + COMMA double* + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, polyval_cmplx, complex_t* + COMMA int + COMMA complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_bartlett, double* x + COMMA int n + COMMA int nfft + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_bartlett_cmplx, complex_t* x + COMMA int n + COMMA int nfft + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_periodogram, double* x + COMMA int n + COMMA int win_type + COMMA double win_param + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_periodogram_cmplx, complex_t* x + COMMA int n + COMMA int win_type + COMMA double win_param + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_welch, double* x + COMMA int n + COMMA int win_type + COMMA double win_param + COMMA int npsd + COMMA int noverlap + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_welch_cmplx, complex_t* x + COMMA int n + COMMA int win_type + COMMA double win_param + COMMA int npsd + COMMA int noverlap + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, randb, double* x + COMMA int n + COMMA random_t* prnd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, randb2, double* x + COMMA int n + COMMA random_t* prnd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, randi, int* x + COMMA int n + COMMA int start + COMMA int stop + COMMA random_t* prnd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, randn, double* x + COMMA int n + COMMA double mu + COMMA double sigma + COMMA random_t* prnd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, randn_cmplx, complex_t* x + COMMA int n + COMMA complex_t* mu + COMMA double sigma + COMMA random_t* prnd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, random_init, random_t* prnd + COMMA int type + COMMA void* seed); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, randu, double* + COMMA int + COMMA random_t* prnd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, ratcompos, double* b + COMMA double* a + COMMA int n + COMMA double* c + COMMA double* d + COMMA int p + COMMA double* beta + COMMA double* alpha); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, re2cmplx, double* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, readbin, char* fn + COMMA void** x + COMMA int* k + COMMA int* dtype); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, signal_pimp, double* + COMMA size_t + COMMA double + COMMA double + COMMA double + COMMA double + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, signal_saw, double* t + COMMA size_t n + COMMA double amp + COMMA double dt + COMMA double period + COMMA double* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, sin_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, sinc, double* x + COMMA int n + COMMA double a + COMMA double* y); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, sine_int, double* x + COMMA int n + COMMA double* si); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, sqrt_cmplx, complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, std, double* x + COMMA int n + COMMA double* s); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, std_cmplx, complex_t* x + COMMA int n + COMMA double* s); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, trapint, double* + COMMA double* + COMMA int + COMMA double*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, trapint_cmplx, double* + COMMA complex_t* + COMMA int + COMMA complex_t*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, unwrap, double* + COMMA int + COMMA double + COMMA double); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, vector_dot, double* x + COMMA double* y + COMMA int n + COMMA double* p); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, verif, double* x + COMMA double* y + COMMA size_t n + COMMA double eps + COMMA double* err); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, verif_data_gen, int len + COMMA int type + COMMA char* fn); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, verif_cmplx, complex_t* x + COMMA complex_t* y + COMMA size_t n + COMMA double eps + COMMA double* err); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(void, verif_str, double* yout + COMMA int nout + COMMA char* str_msg + COMMA char* outfn + COMMA char* logfn); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(void, verif_str_cmplx, complex_t* yout + COMMA int nout + COMMA char* str_msg + COMMA char* outfn + COMMA char* logfn); +/*----------------------------------------------------------------------------*/ +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 + COMMA char*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, writetxt, double* + COMMA double* + COMMA int + COMMA char* ); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, writetxt_3d, double* x + COMMA int nx + COMMA double* y + COMMA int ny + COMMA double* z + COMMA char* fn); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, writetxt_3dline, double* x + COMMA double* y + COMMA double* z + COMMA int n + COMMA char* fn); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, writetxt_cmplx, complex_t* x + COMMA int n + COMMA char* fn); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, writetxt_cmplx_im, double* + COMMA complex_t* + COMMA int + COMMA char*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, writetxt_cmplx_re, double* + COMMA complex_t* + COMMA int + COMMA char*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, writetxt_int, int* + COMMA int* + COMMA int + COMMA char*); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, xcorr, double* x + COMMA int nx + COMMA double* y + COMMA int ny + COMMA int flag + COMMA int nr + COMMA double* r + COMMA double* t); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, xcorr_cmplx, complex_t* x + COMMA int nx + COMMA complex_t* y + COMMA int ny + COMMA int flag + COMMA int nr + COMMA complex_t* r + COMMA double* t); +/*----------------------------------------------------------------------------*/ + + +#ifdef __cplusplus + } +#endif + + +void* dspl_load(); +void dspl_free(void* handle); + + + +#endif /* DSPL_H */ + diff --git a/dspl/dox/en/mainpage.dox b/dspl/dox/en/mainpage.dox index 7da482e..ae34694 100644 --- a/dspl/dox/en/mainpage.dox +++ b/dspl/dox/en/mainpage.dox @@ -57,6 +57,7 @@ Dsplib toolchain includes GCC, GNUPLOT, CodeBlocks IDE, Far file manager and als Digital spectral analysis: \n - \ref DFT_GROUP \n - \ref WIN_GROUP \n + - \ref PSD_GROUP \n - \ref HILBERT_GROUP \n \n diff --git a/dspl/dox/ru/mainpage.dox b/dspl/dox/ru/mainpage.dox index 65bc52b..3f26090 100644 --- a/dspl/dox/ru/mainpage.dox +++ b/dspl/dox/ru/mainpage.dox @@ -49,6 +49,7 @@ DSPL-2.0 --- свободная библиотека алгоритмов циф - \ref DFT_GROUP \n - \ref WIN_GROUP \n + - \ref PSD_GROUP \n - \ref HILBERT_GROUP \n \n diff --git a/dspl/src/psd.c b/dspl/src/psd.c index eb05426..b80b93a 100644 --- a/dspl/src/psd.c +++ b/dspl/src/psd.c @@ -134,33 +134,139 @@ exit_label: + + +#ifdef DOXYGEN_ENGLISH + +#endif +#ifdef DOXYGEN_RUSSIAN + +#endif +int DSPL_API psd_bartlett_cmplx(complex_t* x, int n, int nfft, + fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) +{ + int err, pos, k; + double *pdgr = NULL; + complex_t *tmp = NULL; + fft_t *ptr_fft = NULL; + + pos = 0; + + pdgr = (double*)malloc(nfft * sizeof(double)); + if(!pdgr) + return ERROR_MALLOC; + + if(!pfft) + { + ptr_fft = (fft_t*)malloc(sizeof(fft_t)); + memset(ptr_fft, 0, sizeof(fft_t)); + } + else + ptr_fft = pfft; + + memset(ppsd, 0, nfft * sizeof(double)); + while(pos + nfft <= n) + { + err = fft_mag_cmplx(x + pos, nfft, ptr_fft, fs, + flag & DSPL_FLAG_FFT_SHIFT, pdgr, NULL); + if(err != RES_OK) + goto exit_label; + for(k = 0; k < nfft; k++) + ppsd[k] += pdgr[k]; + pos += nfft; + } + + if(pos < n) + { + tmp = (complex_t*)malloc(nfft * sizeof(complex_t)); + if(!tmp) + { + err = ERROR_MALLOC; + goto exit_label; + } + memset(tmp ,0, nfft * sizeof(complex_t)); + memcpy(tmp, x + pos, (n - pos)*sizeof(complex_t)); + + err = fft_mag_cmplx(tmp, nfft, ptr_fft, fs, + flag & DSPL_FLAG_FFT_SHIFT, pdgr, NULL); + if(err != RES_OK) + goto exit_label; + + for(k = 0; k < nfft; k++) + ppsd[k] += pdgr[k]; + } + + /* fill frequency */ + if(pfrq) + { + if(flag & DSPL_FLAG_FFT_SHIFT) + if(n%2) + err = linspace(-fs*0.5 + fs*0.5/(double)nfft, + fs*0.5 - fs*0.5/(double)nfft, + n, DSPL_SYMMETRIC, pfrq); + else + err = linspace(-fs*0.5, fs*0.5, nfft, DSPL_PERIODIC, pfrq); + else + err = linspace(0, fs, nfft, DSPL_PERIODIC, pfrq); + } + + /* scale magnitude */ + if(flag & DSPL_FLAG_LOGMAG) + { + for(k = 0; k < nfft; k++) + ppsd[k] = 10.0 * log10(ppsd[k] / (double)n / fs); + } + else + { + for(k = 0; k < nfft; k++) + ppsd[k] /= (double)n * fs; + } + + +exit_label: + if(pdgr) + free(pdgr); + if(tmp) + free(tmp); + if(ptr_fft && (ptr_fft != pfft)) + { + fft_free(ptr_fft); + free(ptr_fft); + } + return err; +} + + + + + + #ifdef DOXYGEN_ENGLISH #endif #ifdef DOXYGEN_RUSSIAN /*! **************************************************************************** \ingroup PSD_GROUP -\fn -int psd_periodogram(double* x, int n, - int win_type, double win_param, - fft_t* pfft, double fs, - int flag, double* ppsd, double* pfrq) -\brief Непараметрическая оценка спектральной плотности сигнала методом - модифицированной периодограммы +\fn int psd_periodogram(double* x, int n, + int win_type, double win_param, + fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) +\brief Непараметрическая оценка спектральной плотности мощности (СПМ) +вещественного сигнала методом модифицированной периодограммы. Функция рассчитывает спектральную плотность мощности \f$ X(f) \f$ выборки сигнала длительности \$n \$ отсчетов методом модифицированной периодограммы: \f[ - X(f) = \frac{1}{N U F_s} \Big| \sum_{m = 0}^{n-1} w(m) x(m) \exp - \left( -j \frac{2\pi} f m \right) \Big|^2, + X(f) = \frac{1}{U F_s} \left| \sum_{m = 0}^{n-1} w(m) x(m) \exp + \left( -j 2\pi f m \right) \right|^2, \f] где \f$ w(m) \f$ -- отсчёты оконной функции, \f$ F_s \f$ -- частота дискретизации (Гц), \f$ U \f$ нормировочный коэффициент равный - \f[ - U \sum_{m = 0}^{n-1} w^2(m) + U = \sum_{m = 0}^{n-1} w^2(m) \f] При использовании прямоугольного окна модифицированная периодограмма переходит @@ -171,12 +277,9 @@ int psd_periodogram(double* x, int n, (по умолчанию), или от \f$-F_s /2 \f$ до \f$F_s /2 \f$, если установлен флаг расчета двусторонней периодограммы. - -\note Периодограмма, как стандартная, так и модифицированная возвращает -асимптотически несмещенную, но несостоятельную оценку СПМ (уроверь флуктуаций -шумовой составляющей СПМ не уменьшается с ростом длины выборки `n`). - - +\note Периодограмма возвращает асимптотически несмещенную, +но несостоятельную оценку СПМ (уровень флуктуаций шумовой составляющей СПМ +не уменьшается с ростом длины выборки `n`). \param[in] x Указатель на входной вектор вещественного сигнала \f$x(m)\f$, @@ -200,11 +303,12 @@ int psd_periodogram(double* x, int n, Подробнее смотри описание функции \ref window. \n\n \param[in] pfft -Указатель на структуру `fft_t`. \n +Указатель на структуру \ref fft_t. \n Указатель может быть `NULL`. В этом случае объект структуры будет создан внутри функции и удален перед завершением.\n Если предполагается многократный вызов функции, то рекомендуется создать -объект `fft_t` и передавать в функцию, чтобы не создавать его каждый раз. \n\n +объект \ref fft_t и передавать в функцию, чтобы не +создавать его каждый раз. \n\n \param[in] fs частота дискретизации выборки исходного сигнала (Гц). \n\n @@ -217,25 +321,50 @@ DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2) \endverbatim \param[in, out] ppsd -Указатель на вектор СПМ расчитанных по входному сигналу $x$. \n +Указатель на вектор СПМ рассчитанных по входному сигналу $x$. \n Размер вектора `[n x 1]`. \n Память должна быть выделена. \n\n \param[in, out] pfrq -Указатель на вектор частоты, соответсвующей +Указатель на вектор частоты, соответствующей значениям рассчитанного вектора СПМ. \n Размер вектора `[n x 1]`. \n -Указатель можеть быть `NULL`,в этом случае вектор частоты не +Указатель может быть `NULL`,в этом случае вектор частоты не рассчитывается и не возвращается. \n\n - \return `RES_OK` если расчет произведен успешно. \n В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n \n +Пример периодограммных оценок СПМ для различной длины выборки сигнала: + +\include psd_periodogram_test.c + +Программа производит расчет СПМ сигнала, состоящего из двух гармоник на +фоне белого гауссова шума. Расчет ведется по выборкам длины 128, 1024 и +8192 отсчетов. + +В результате периодограммы (стандартная с прямоугольным окном +и модифицированная с окном Блэкмана) выводятся на графики: + +`n = 8192` точек (черная -- классическая периодограмма с прямоугольным окном, +зеленая -- модифицированная с окном Блэкмана): +\image html psd_perodogram_8192.png + +`n = 1024` точек (черная -- классическая периодограмма с прямоугольным окном, +зеленая -- модифицированная с окном Блэкмана): +\image html psd_perodogram_1024.png + +`n = 128` точек (черная -- классическая периодограмма с прямоугольным окном, +зеленая -- модифицированная с окном Блэкмана): +\image html psd_perodogram_128.png + +Можно видеть, что модифицированная периодограмма позволяет снизить +растекание СПМ, однако уровень флуктуация шума не уменьшается с увеличением +размера выборки от 128 до 8192 отсчетов (оценка несостоятельная). + \author Бахурин Сергей www.dsplib.org ***************************************************************************** */ - #endif int DSPL_API psd_periodogram(double* x, int n, int win_type, double win_param, @@ -347,7 +476,124 @@ exit_label: #endif #ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup PSD_GROUP +\fn int psd_periodogram_cmplx(complex_t* x, int n, + int win_type, double win_param, + fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) +\brief Непараметрическая оценка спектральной плотности мощности (СПМ) +комплексного сигнала методом модифицированной периодограммы. +Функция рассчитывает спектральную плотность мощности \f$ X(f) \f$ +выборки сигнала длительности \$n \$ отсчетов методом модифицированной +периодограммы: + +\f[ + X(f) = \frac{1}{U F_s} \left| \sum_{m = 0}^{n-1} w(m) x(m) \exp + \left( -j 2\pi f m \right) \right|^2, +\f] +где \f$ w(m) \f$ -- отсчёты оконной функции, \f$ F_s \f$ -- частота +дискретизации (Гц), \f$ U \f$ нормировочный коэффициент равный +\f[ + U = \sum_{m = 0}^{n-1} w^2(m) +\f] + +При использовании прямоугольного окна модифицированная периодограмма переходит +в стандартную периодограмму. + +Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого +преобразования Фурье, для дискретной сетки частот от 0 Гц до \f$ F_s \f$ Гц +(по умолчанию), или от \f$-F_s /2 \f$ до \f$F_s /2 \f$, если установлен флаг +расчета двусторонней периодограммы. + +\note Периодограмма возвращает асимптотически несмещенную, +но несостоятельную оценку СПМ (уровень флуктуаций шумовой составляющей СПМ +не уменьшается с ростом длины выборки `n`). + +\param[in] x +Указатель на входной вектор вещественного сигнала \f$x(m)\f$, +\f$ m = 0 \ldots n-1 \f$. \n +Размер вектора `[n x 1]`. \n \n + +\param[in] n +Размер вектора входного сигнала. +Также размер выходного вектора СПМ и +вектора частоты также равны `n`.\n\n + +\param[in] win_type +Тип оконной функции, применяемой для модифицированной периодограммы.\n +Подробнее смотри описание функции \ref window. \n\n + + +\param[in] win_type +Параметр оконной функции.\n +Данный параметр используется, если задано параметрическая оконная функция. +Для непараметрических окон данный параметр игнорируется.\n +Подробнее смотри описание функции \ref window. \n\n + +\param[in] pfft +Указатель на структуру \ref fft_t. \n +Указатель может быть `NULL`. В этом случае объект структуры будет +создан внутри функции и удален перед завершением.\n +Если предполагается многократный вызов функции, то рекомендуется создать +объект \ref fft_t и передавать в функцию, чтобы не создавать его каждый раз. \n\n + +\param[in] fs +частота дискретизации выборки исходного сигнала (Гц). \n\n + +\param[in] flag +Комбинация битовых флагов, задающих режим расчета: +\verbatim +DSPL_FLAG_LOGMAG - СПМ считать в логарифмическом масштабе в единицах дБ/Гц +DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2) +\endverbatim + +\param[in, out] ppsd +Указатель на вектор СПМ рассчитанных по входному сигналу $x$. \n +Размер вектора `[n x 1]`. \n +Память должна быть выделена. \n\n + +\param[in, out] pfrq +Указатель на вектор частоты, соответствующей +значениям рассчитанного вектора СПМ. \n +Размер вектора `[n x 1]`. \n +Указатель может быть `NULL`,в этом случае вектор частоты не +рассчитывается и не возвращается. \n\n + +\return +`RES_OK` если расчет произведен успешно. \n +В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n \n + +Пример периодограммных оценок СПМ для различной длины выборки сигнала: + +\include psd_periodogram_cmplx_test.c + +Программа производит расчет СПМ сигнала, состоящего из двух гармоник на +фоне белого гауссова шума. Расчет ведется по выборкам длины 128, 1024 и +8192 отсчетов. + +В результате периодограммы (стандартная с прямоугольным окном +и модифицированная с окном Блэкмана) выводятся на графики: + +`n = 8192` точек (черная -- классическая периодограмма с прямоугольным окном, +зеленая -- модифицированная с окном Блэкмана): +\image html psd_perodogram_cmplx_8192.png + +`n = 1024` точек (черная -- классическая периодограмма с прямоугольным окном, +зеленая -- модифицированная с окном Блэкмана): +\image html psd_perodogram_cmplx_1024.png + +`n = 128` точек (черная -- классическая периодограмма с прямоугольным окном, +зеленая -- модифицированная с окном Блэкмана): +\image html psd_perodogram_cmplx_128.png + +Можно видеть, что модифицированная периодограмма позволяет снизить +растекание СПМ, однако уровень флуктуация шума не уменьшается с увеличением +размера выборки от 128 до 8192 отсчетов (оценка несостоятельная). + +\author Бахурин Сергей www.dsplib.org +***************************************************************************** */ #endif int DSPL_API psd_periodogram_cmplx(complex_t* x, int n, int win_type, double win_param, diff --git a/dspl/src/randgen.c b/dspl/src/randgen.c index fc0be1d..2d987a8 100644 --- a/dspl/src/randgen.c +++ b/dspl/src/randgen.c @@ -627,6 +627,35 @@ exit_label: +#ifdef DOXYGEN_ENGLISH + +#endif +#ifdef DOXYGEN_RUSSIAN + +#endif +int DSPL_API randn_cmplx(complex_t* x, int n, complex_t* mu, + double sigma, random_t* prnd) +{ + int err, i; + + err = randn((double*)x, 2*n, 0.0, sigma / M_SQRT2, prnd); + if(err!= RES_OK) + return err; + if(mu) + { + for(i = 0; i < n; i++) + { + RE(x[i]) += RE(mu[0]); + IM(x[i]) += IM(mu[0]); + } + } + return err; +} + + + + + #ifdef DOXYGEN_ENGLISH diff --git a/examples/src/psd_bartlett_cmplx_test.c b/examples/src/psd_bartlett_cmplx_test.c new file mode 100644 index 0000000..fa989bf --- /dev/null +++ b/examples/src/psd_bartlett_cmplx_test.c @@ -0,0 +1,95 @@ +#include +#include +#include +#include "dspl.h" + + +#define N 8192 +#define FS 1.0 + + + +void psd_plot(int argc, char* argv[], char* fn_png, char* fn_data) +{ + char cmd[512] = {0}; + void* hplot; /* GNUPLOT handle */ + /* plotting by GNUPLOT */ + gnuplot_create(argc, argv, 680, 480, fn_png, &hplot); + gnuplot_cmd(hplot, "unset key"); + gnuplot_cmd(hplot, "set grid"); + gnuplot_cmd(hplot, "set xlabel 'frequency, Hz'"); + gnuplot_cmd(hplot, "set ylabel 'Power Spectral Density, [dB/Hz]'"); + gnuplot_cmd(hplot, "set yrange [-60: 40]"); + + sprintf(cmd, "plot '%s' w l lt -1", fn_data); + + gnuplot_cmd(hplot, cmd); + gnuplot_close(hplot); +} + + +int main(int argc, char* argv[]) +{ + void* hdspl; /* DSPL handle */ + + hdspl = dspl_load(); /* Load DSPL function */ + random_t rnd = {0}; /* random structure */ + + complex_t *x = NULL; + double *psd = NULL; + double *frq = NULL; + + int k, err; + + x = (complex_t*)malloc(N*sizeof(complex_t)); + psd = (double*)malloc(N*sizeof(double)); + frq = (double*)malloc(N*sizeof(double)); + + /* input signal fill as noise -20 dB/Hz power spectrum density */ + random_init(&rnd, RAND_TYPE_MRG32K3A, NULL); + randn_cmplx(x, N, NULL, 0.1, &rnd); + + /* x[k] = 0.1 * cos(2 * pi * 0.2 * k) + cos(2 * pi * 0.26 * k) + noise */ + for(k = 0; k < N; k++) + { + RE(x[k]) += cos(M_2PI * 0.26 * (double)k) + + 0.1*cos(M_2PI*0.2* (double)k); + + IM(x[k]) += sin(M_2PI * 0.26 * (double)k) + + 0.1*sin(M_2PI*0.2* (double)k); + } + /* Twosided PSD logscale magnitude n = 8192, nfft = 8192. + This case is periodogram */ + err = psd_bartlett_cmplx(x, N, N, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N, "dat/psd_bartlett_cmplx_8192.txt"); + psd_plot(argc, argv, "img/psd_bartlett_cmplx_8192.png", + "dat/psd_bartlett_cmplx_8192.txt"); + + + /* Twosided PSD with logscale magnitude n = 8192, nfft = 1024 */ + err = psd_bartlett_cmplx(x, N, N/8, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/8, "dat/psd_bartlett_cmplx_1024.txt"); + psd_plot(argc, argv, "img/psd_bartlett_cmplx_1024.png", + "dat/psd_bartlett_cmplx_1024.txt"); + + + /* Twosided PSD with logscale magnitude n = 8192, nfft = 128 */ + err = psd_bartlett_cmplx(x, N, N/64, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/64, "dat/psd_bartlett_cmplx_128.txt"); + psd_plot(argc, argv, "img/psd_bartlett_cmplx_128.png", + "dat/psd_bartlett_cmplx_128.txt"); + + dspl_free(hdspl); /* free dspl handle */ + if(x) + free(x); + if(psd) + free(psd); + if(frq) + free(frq); + return 0; +} + + diff --git a/examples/src/psd_bartlett_test.c b/examples/src/psd_bartlett_test.c new file mode 100644 index 0000000..ab49576 --- /dev/null +++ b/examples/src/psd_bartlett_test.c @@ -0,0 +1,90 @@ +#include +#include +#include +#include "dspl.h" + + +#define N 8192 +#define FS 1.0 + + + +void psd_plot(int argc, char* argv[], char* fn_png, char* fn_data) +{ + char cmd[512] = {0}; + void* hplot; /* GNUPLOT handle */ + /* plotting by GNUPLOT */ + gnuplot_create(argc, argv, 680, 480, fn_png, &hplot); + gnuplot_cmd(hplot, "unset key"); + gnuplot_cmd(hplot, "set grid"); + gnuplot_cmd(hplot, "set xlabel 'frequency, Hz'"); + gnuplot_cmd(hplot, "set ylabel 'Power Spectral Density, [dB/Hz]'"); + gnuplot_cmd(hplot, "set yrange [-60: 40]"); + + sprintf(cmd, "plot '%s' w l lt -1", fn_data); + + gnuplot_cmd(hplot, cmd); + gnuplot_close(hplot); +} + + +int main(int argc, char* argv[]) +{ + void* hdspl; /* DSPL handle */ + + hdspl = dspl_load(); /* Load DSPL function */ + random_t rnd = {0}; /* random structure */ + + double *x = NULL; + double *psd = NULL; + double *frq = NULL; + + int k, err; + + x = (double*)malloc(N*sizeof(double)); + psd = (double*)malloc(N*sizeof(double)); + frq = (double*)malloc(N*sizeof(double)); + + /* input signal fill as noise -20 dB/Hz power spectrum density */ + random_init(&rnd, RAND_TYPE_MRG32K3A, NULL); + randn(x, N, 0.0, 0.1, &rnd); + + /* x[k] = 0.1 * cos(2 * pi * 0.2 * k) + cos(2 * pi * 0.26 * k) + noise */ + for(k = 0; k < N; k++) + x[k] += cos(M_2PI * 0.26 * (double)k) + 0.1*cos(M_2PI*0.2* (double)k); + + /* Twosided PSD logscale magnitude n = 8192, nfft = 8192. + This case is periodogram */ + err = psd_bartlett(x, N, N, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N, "dat/psd_bartlett_8192.txt"); + psd_plot(argc, argv, "img/psd_bartlett_8192.png", + "dat/psd_bartlett_8192.txt"); + + + /* Twosided PSD with logscale magnitude n = 8192, nfft = 1024 */ + err = psd_bartlett(x, N, N/8, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/8, "dat/psd_bartlett_1024.txt"); + psd_plot(argc, argv, "img/psd_bartlett_1024.png", + "dat/psd_bartlett_1024.txt"); + + + /* Twosided PSD with logscale magnitude n = 8192, nfft = 128 */ + err = psd_bartlett(x, N, N/64, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/64, "dat/psd_bartlett_128.txt"); + psd_plot(argc, argv, "img/psd_bartlett_128.png", + "dat/psd_bartlett_128.txt"); + + dspl_free(hdspl); /* free dspl handle */ + if(x) + free(x); + if(psd) + free(psd); + if(frq) + free(frq); + return 0; +} + + diff --git a/examples/src/psd_periodogram_cmplx_test.c b/examples/src/psd_periodogram_cmplx_test.c new file mode 100644 index 0000000..b6bddb1 --- /dev/null +++ b/examples/src/psd_periodogram_cmplx_test.c @@ -0,0 +1,116 @@ +#include +#include +#include +#include "dspl.h" + + +#define N 8192 +#define FS 1.0 + +void psd_plot(int argc, char* argv[], char* fn_png, + char* fn_data0, char* fn_data1) +{ + char cmd[512] = {0}; + void* hplot; /* GNUPLOT handle */ + /* plotting by GNUPLOT */ + gnuplot_create(argc, argv, 680, 480, fn_png, &hplot); + gnuplot_cmd(hplot, "unset key"); + gnuplot_cmd(hplot, "set grid"); + gnuplot_cmd(hplot, "set xlabel 'frequency, Hz'"); + gnuplot_cmd(hplot, "set ylabel 'Power Spectral Density, [dB/Hz]'"); + gnuplot_cmd(hplot, "set yrange [-60: 40]"); + + sprintf(cmd, "plot '%s' w l lt -1,\\", fn_data0); + gnuplot_cmd(hplot, cmd); + + memset(cmd, 0, 512); + sprintf(cmd, "'%s' w l lt 2", fn_data1); + gnuplot_cmd(hplot, cmd); + + gnuplot_close(hplot); +} + + +int main(int argc, char* argv[]) +{ + void* hdspl; /* DSPL handle */ + void* hplot; /* GNUPLOT handle */ + hdspl = dspl_load(); // Load DSPL function + random_t rnd = {0}; /* random structure */ + + complex_t *x = NULL; + double *psd = NULL; + double *frq = NULL; + + int k, err; + + x = (complex_t*)malloc(N*sizeof(complex_t)); + psd = (double*)malloc(N*sizeof(double)); + frq = (double*)malloc(N*sizeof(double)); + + /* input signal fill as noise -40 dB/Hz power spectrum density */ + random_init(&rnd, RAND_TYPE_MRG32K3A, NULL); + randn_cmplx(x, N, NULL, 0.1, &rnd); + + /* x[k] = 0.1*cos(2*pi*0.1*k) + cos(2*pi*0.2*k) + noise */ + for(k = 0; k < N; k++) + { + RE(x[k]) += cos(M_2PI * 0.26 * (double)k) + + 0.1*cos(M_2PI * 0.20 * (double)k); + + IM(x[k]) += sin(M_2PI * 0.26 * (double)k) + + 0.1*sin(M_2PI * 0.20 * (double)k); + } + /* Twosided PSD with logscale magnitude 8192 points with rect window */ + err = psd_periodogram_cmplx(x, N, DSPL_WIN_RECT, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N, "dat/psd_periodogram_8192_win_rect.txt"); + + /* Twosided PSD with logscale magnitude 8192 points with blackman window */ + err = psd_periodogram_cmplx(x, N, DSPL_WIN_BLACKMAN, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N, "dat/psd_periodogram_8192_win_blackman.txt"); + + psd_plot(argc, argv, "img/psd_perodogram_cmplx_8192.png", + "dat/psd_periodogram_8192_win_rect.txt", + "dat/psd_periodogram_8192_win_blackman.txt"); + + /* Twosided PSD with logscale magnitude 1024 points with rect window */ + err = psd_periodogram_cmplx(x, N/8, DSPL_WIN_RECT, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/8, "dat/psd_periodogram_1024_win_rect.txt"); + + /* Twosided PSD with logscale magnitude 1024 points with blackman window */ + err = psd_periodogram_cmplx(x, N/8, DSPL_WIN_BLACKMAN, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/8, "dat/psd_periodogram_1024_win_blackman.txt"); + + psd_plot(argc, argv, "img/psd_perodogram_cmplx_1024.png", + "dat/psd_periodogram_1024_win_rect.txt", + "dat/psd_periodogram_1024_win_blackman.txt"); + + /* Twosided PSD with logscale magnitude 128 points with rect window */ + err = psd_periodogram_cmplx(x, N/64, DSPL_WIN_RECT, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/64, "dat/psd_periodogram_128_win_rect.txt"); + + /* Twosided PSD with logscale magnitude 128 points with blackman window */ + err = psd_periodogram_cmplx(x, N/64, DSPL_WIN_BLACKMAN, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/64, "dat/psd_periodogram_128_win_blackman.txt"); + + psd_plot(argc, argv, "img/psd_perodogram_cmplx_128.png", + "dat/psd_periodogram_128_win_rect.txt", + "dat/psd_periodogram_128_win_blackman.txt"); + + dspl_free(hdspl); // free dspl handle + if(x) + free(x); + if(psd) + free(psd); + if(frq) + free(frq); + return 0; +} + + diff --git a/examples/src/psd_periodogram_test.c b/examples/src/psd_periodogram_test.c new file mode 100644 index 0000000..a54cd6f --- /dev/null +++ b/examples/src/psd_periodogram_test.c @@ -0,0 +1,119 @@ +#include +#include +#include +#include "dspl.h" + + +#define N 8192 +#define FS 1.0 + + +void psd_plot(int argc, char* argv[], char* fn_png, + char* fn_data0, char* fn_data1) +{ + char cmd[512] = {0}; + void* hplot; /* GNUPLOT handle */ + /* plotting by GNUPLOT */ + gnuplot_create(argc, argv, 680, 480, fn_png, &hplot); + gnuplot_cmd(hplot, "unset key"); + gnuplot_cmd(hplot, "set grid"); + gnuplot_cmd(hplot, "set xlabel 'frequency, Hz'"); + gnuplot_cmd(hplot, "set ylabel 'Power Spectral Density, [dB/Hz]'"); + gnuplot_cmd(hplot, "set yrange [-60: 40]"); + + sprintf(cmd, "plot '%s' w l lt -1,\\", fn_data0); + gnuplot_cmd(hplot, cmd); + + memset(cmd, 0, 512); + sprintf(cmd, "'%s' w l lt 2", fn_data1); + gnuplot_cmd(hplot, cmd); + + gnuplot_close(hplot); +} + + +int main(int argc, char* argv[]) +{ + void* hdspl; /* DSPL handle */ + void* hplot; /* GNUPLOT handle */ + hdspl = dspl_load(); // Load DSPL function + random_t rnd = {0}; /* random structure */ + + double *x = NULL; + double *psd = NULL; + double *frq = NULL; + + int k, err; + + x = (double*)malloc(N*sizeof(double)); + psd = (double*)malloc(N*sizeof(double)); + frq = (double*)malloc(N*sizeof(double)); + + /* input signal fill as noise -40 dB/Hz power spectrum density */ + random_init(&rnd, RAND_TYPE_MRG32K3A, NULL); + randn(x, N, 0.0, 0.1, &rnd); + + /* x[k] = 0.1*cos(2*pi*0.1*k) + cos(2*pi*0.2*k) + noise */ + for(k = 0; k < N; k++) + x[k] += cos(M_2PI * 0.26 * (double)k) + 0.1*cos(M_2PI*0.2* (double)k); + + + + /* Twosided PSD with logscale magnitude 8192 points with rect window */ + err = psd_periodogram(x, N, DSPL_WIN_RECT, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N, "dat/psd_periodogram_8192_win_rect.txt"); + + + /* Twosided PSD with logscale magnitude 8192 points with blackman window */ + err = psd_periodogram(x, N, DSPL_WIN_BLACKMAN, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N, "dat/psd_periodogram_8192_win_blackman.txt"); + + psd_plot(argc, argv, "img/psd_perodogram_8192.png", + "dat/psd_periodogram_8192_win_rect.txt", + "dat/psd_periodogram_8192_win_blackman.txt"); + + + /* Twosided PSD with logscale magnitude 1024 points with rect window */ + err = psd_periodogram(x, N/8, DSPL_WIN_RECT, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/8, "dat/psd_periodogram_1024_win_rect.txt"); + + /* Twosided PSD with logscale magnitude 1024 points with blackman window */ + err = psd_periodogram(x, N/8, DSPL_WIN_BLACKMAN, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/8, "dat/psd_periodogram_1024_win_blackman.txt"); + + psd_plot(argc, argv, "img/psd_perodogram_1024.png", + "dat/psd_periodogram_1024_win_rect.txt", + "dat/psd_periodogram_1024_win_blackman.txt"); + + + /* Twosided PSD with logscale magnitude 128 points with rect window */ + err = psd_periodogram(x, N/64, DSPL_WIN_RECT, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/64, "dat/psd_periodogram_128_win_rect.txt"); + + /* Twosided PSD with logscale magnitude 128 points with blackman window */ + err = psd_periodogram(x, N/64, DSPL_WIN_BLACKMAN, 0, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, N/64, "dat/psd_periodogram_128_win_blackman.txt"); + + psd_plot(argc, argv, "img/psd_perodogram_128.png", + "dat/psd_periodogram_128_win_rect.txt", + "dat/psd_periodogram_128_win_blackman.txt"); + + + + dspl_free(hdspl); // free dspl handle + if(x) + free(x); + if(psd) + free(psd); + if(frq) + free(frq); + return 0; +} + + diff --git a/include/dspl.c b/include/dspl.c index 29b9f8d..4846af8 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -154,7 +154,9 @@ p_polyroots polyroots ; p_polyval polyval ; p_polyval_cmplx polyval_cmplx ; p_psd_bartlett psd_bartlett ; +p_psd_bartlett_cmplx psd_bartlett_cmplx ; p_psd_periodogram psd_periodogram ; +p_psd_periodogram_cmplx psd_periodogram_cmplx ; p_psd_welch psd_welch ; p_psd_welch_cmplx psd_welch_cmplx ; @@ -162,6 +164,7 @@ p_randb randb ; p_randb2 randb2 ; p_randi randi ; p_randn randn ; +p_randn_cmplx randn_cmplx ; p_random_init random_init ; p_randu randu ; p_ratcompos ratcompos ; @@ -472,7 +475,9 @@ void* dspl_load() LOAD_FUNC(polyval); LOAD_FUNC(polyval_cmplx); LOAD_FUNC(psd_bartlett); + LOAD_FUNC(psd_bartlett_cmplx); LOAD_FUNC(psd_periodogram); + LOAD_FUNC(psd_periodogram_cmplx); LOAD_FUNC(psd_welch); LOAD_FUNC(psd_welch_cmplx); @@ -480,6 +485,7 @@ void* dspl_load() LOAD_FUNC(randb); LOAD_FUNC(randb2); LOAD_FUNC(randn); + LOAD_FUNC(randn_cmplx); LOAD_FUNC(random_init); LOAD_FUNC(randu); LOAD_FUNC(ratcompos); diff --git a/include/dspl.h b/include/dspl.h index 507eb25..b5b7bef 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -1298,6 +1298,15 @@ DECLARE_FUNC(int, psd_bartlett, double* x COMMA double* ppsd COMMA double* pfrq); /*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_bartlett_cmplx, complex_t* x + COMMA int n + COMMA int nfft + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ DECLARE_FUNC(int, psd_periodogram, double* x COMMA int n COMMA int win_type @@ -1308,6 +1317,16 @@ DECLARE_FUNC(int, psd_periodogram, double* x COMMA double* ppsd COMMA double* pfrq); /*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, psd_periodogram_cmplx, complex_t* x + COMMA int n + COMMA int win_type + COMMA double win_param + COMMA fft_t* pfft + COMMA double fs + COMMA int flag + COMMA double* ppsd + COMMA double* pfrq); +/*----------------------------------------------------------------------------*/ DECLARE_FUNC(int, psd_welch, double* x COMMA int n COMMA int win_type @@ -1346,10 +1365,16 @@ DECLARE_FUNC(int, randi, int* x COMMA int stop COMMA random_t* prnd); /*----------------------------------------------------------------------------*/ -DECLARE_FUNC(int, randn, double* - COMMA int - COMMA double - COMMA double +DECLARE_FUNC(int, randn, double* x + COMMA int n + COMMA double mu + COMMA double sigma + COMMA random_t* prnd); +/*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, randn_cmplx, complex_t* x + COMMA int n + COMMA complex_t* mu + COMMA double sigma COMMA random_t* prnd); /*----------------------------------------------------------------------------*/ DECLARE_FUNC(int, random_init, random_t* prnd