diff --git a/dspl/dox/ru/mainpage.dox b/dspl/dox/ru/mainpage.dox index 4fb38e1..65bc52b 100644 --- a/dspl/dox/ru/mainpage.dox +++ b/dspl/dox/ru/mainpage.dox @@ -14,7 +14,7 @@ DSPL-2.0 --- свободная библиотека алгоритмов циф -Исходные коды библиотеки libdspl-2.0 на GitHub +libdspl-2.0 на GitHub diff --git a/dspl/src/fft.c b/dspl/src/fft.c index a9af427..259299c 100644 --- a/dspl/src/fft.c +++ b/dspl/src/fft.c @@ -21,6 +21,7 @@ #include #include #include +#include #include "dspl.h" #include "dspl_internal.h" @@ -915,7 +916,6 @@ int DSPL_API fft_mag(double* x, int n, fft_t* pfft, { int k, err = RES_OK; complex_t *X = NULL; - if(!x || !pfft) return ERROR_PTR; @@ -931,7 +931,7 @@ int DSPL_API fft_mag(double* x, int n, fft_t* pfft, if(flag & DSPL_FLAG_LOGMAG) for(k = 0; k < n; k++) - mag[k] = 10.0*log10(ABSSQR(X[k])); + mag[k] = 10.0*log10(ABSSQR(X[k])+DBL_EPSILON); else for(k = 0; k < n; k++) mag[k] = ABS(X[k]); diff --git a/dspl/src/inout.c b/dspl/src/inout.c index 510362f..7edd086 100644 --- a/dspl/src/inout.c +++ b/dspl/src/inout.c @@ -505,7 +505,10 @@ int DSPL_API writetxt(double* x, double* y, int n, char* fn) if(y) for(k = 0; k < n; k++) - fprintf(pFile, "%+.12E\t%+.12E\n", x[k], y[k]); + if(!isnan(x[k]) && !isnan(y[k]) && !isinf(x[k]) && !isinf(y[k])) + fprintf(pFile, "%+.12E\t%+.12E\n", x[k], y[k]); + else + break; else for(k = 0; k < n; k++) diff --git a/dspl/src/win.c b/dspl/src/win.c index c70bf96..1c2886a 100644 --- a/dspl/src/win.c +++ b/dspl/src/win.c @@ -28,7 +28,141 @@ #ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup WIN_GROUP +\fn int window(double* w, int n, int win_type, double param) +\brief Window function calculation +The function calculates a periodic or symmetric window function +according to parameter `win_type`. \n + +A periodic window function is used for spectral analysis, +and a symmetric window function can be used to design FIR filters. + +\param [in,out] w +Pointer to the window. \n +Vector size is `[n x 1]`. \n +Memory must be allocated. \n +The calculated window function will be placed at the given address. \n \n + +\param [in] n +Size of window function `w` vector. \n \n + +\param [in] win_type +Combination of flags for specifying the type of window function. \n +Combination of `DSPL_WIN_MASK | DSPL_WIN_SYM_MASK` bit masks +is used to set the window type.\n +Bit mask `DSPL_WIN_MASK` sets the window type. +Can be one of follow: \n +\verbatim +------------------------------------------------------------------------- + win_type | Description +-----------------------------|------------------------------------------- + DSPL_WIN_BARTLETT | Nonparametric Bartlett window +-----------------------------|------------------------------------------- + DSPL_WIN_BARTLETT_HANN | Nonparametric Bartlett-Hann window +-----------------------------|------------------------------------------- + DSPL_WIN_BLACKMAN | Nonparametric Blackman window +-----------------------------|------------------------------------------- + DSPL_WIN_BLACKMAN_HARRIS | Nonparametric Blackman-Harris window +-----------------------------|------------------------------------------- + DSPL_WIN_BLACKMAN_NUTTALL | Nonparametric Blackman-Nuttall +-----------------------------|------------------------------------------- + DSPL_WIN_CHEBY | Parametric Dolph-Chebyshev window. + | Parametr `win_param` sets sidelobe attenuation + | level in dB. +-----------------------------|------------------------------------------- + DSPL_WIN_COS | Nonparametric Cosine window +-----------------------------|------------------------------------------- + DSPL_WIN_FLAT_TOP | Nonparametric maxflat window +-----------------------------|------------------------------------------- + DSPL_WIN_GAUSSIAN | Nonparametric Gauss window +-----------------------------|------------------------------------------- + DSPL_WIN_HAMMING | Nonparametric Hamming window +-----------------------------|------------------------------------------- + DSPL_WIN_HANN | Nonparametric Hann window +-----------------------------|------------------------------------------- + DSPL_WIN_KAISER | Parametric Kaiser window +-----------------------------|------------------------------------------- + DSPL_WIN_LANCZOS | Nonparametric Lanczos window +-----------------------------|------------------------------------------- + DSPL_WIN_NUTTALL | Nonparametric Nuttall window +-----------------------------|------------------------------------------- + DSPL_WIN_RECT | Nonparametric rectangular window +------------------------------------------------------------------------- +\endverbatim +\n +Bit mask `DSPL_WIN_SYM_MASK` sets window function symmetry: \n +\verbatim +------------------------------------------------------------------------- + DSPL_WIN_SYM_MASK | Description +-----------------------------|------------------------------------------- + DSPL_WIN_SYMMETRIC | Symmetry window (default value) + DSPL_WIN_PERIODIC | Periodic window +------------------------------------------------------------------------- +\endverbatim + \n \n + +\param [in] param +Window function parameter. \n +This parameter is using only to parametric window functions, +and ignored for nonparametric windows. \n +\n + +\return +`RES_OK` if window function is calculated successfully. \n +Else \ref ERROR_CODE_GROUP "error code". + +The following program calculates 64 samples window functions, +draws their spectrum when using the bin indices of +the discrete Fourier transform along the frequency axis. + +\include windows_test.c + +A personal graph is displayed for each type of window function. + +Rectangular window +\image html win_rect.png +\n +\n + +Nonparametric windows +\image html win_bartlett.png +\image html win_flattop.png +\image html win_bartletthann.png +\image html win_hann.png +\image html win_hamming.png +\image html win_lanczos.png +\image html win_blackman.png +\image html win_blackmanharris.png +\image html win_blackmannuttall.png +\image html win_cos.png +\image html win_nuttall.png +\n +\n + +Parametric Dolph-Chebyshev windows +\image html win_cheby50.png +\image html win_cheby80.png +\image html win_cheby120.png +\n +\n + +Parametric Gaussian windows +\image html win_gaussian0p5.png +\image html win_gaussian0p3.png + +\n +\n +Parametric Kaiser windows +\image html win_kaiser4p0.png +\image html win_kaiser8p0.png +\image html win_kaiser12p0.png +\n +\n + +\author Sergey Bakhurin. www.dsplib.org +***************************************************************************** */ #endif #ifdef DOXYGEN_RUSSIAN /*! **************************************************************************** @@ -120,6 +254,54 @@ `RES_OK` если оконная функция рассчитана успешно. \n В противном случае \ref ERROR_CODE_GROUP "код ошибки". +Следующая программа производит расчет оконных функций длительности 64 отсчета, +строит их спектральную плотность при использовании по оси частот индексы бинов +дискретного преобразования Фурье. + +\include windows_test.c + +Для каждого вида оконной функция выводится персональный график. + +Прямоугольное окно +\image html win_rect.png +\n +\n + +Непраметрические окна +\image html win_bartlett.png +\image html win_flattop.png +\image html win_bartletthann.png +\image html win_hann.png +\image html win_hamming.png +\image html win_lanczos.png +\image html win_blackman.png +\image html win_blackmanharris.png +\image html win_blackmannuttall.png +\image html win_cos.png +\image html win_nuttall.png +\n +\n + +Параметрические окна Дольф-Чебышева +\image html win_cheby50.png +\image html win_cheby80.png +\image html win_cheby120.png +\n +\n + +Параметрические окна Гаусса +\image html win_gaussian0p5.png +\image html win_gaussian0p3.png + +\n +\n +Параметрические окна Кайзера +\image html win_kaiser4p0.png +\image html win_kaiser8p0.png +\image html win_kaiser12p0.png +\n +\n + \author Бахурин Сергей. www.dsplib.org ***************************************************************************** */ #endif @@ -484,7 +666,7 @@ Gaussian window function ******************************************************************************/ int win_gaussian(double *w, int n, int win_type, double alpha) { - double x = 0.0; + double a = 0.0; double y; double sigma; int i; @@ -496,16 +678,16 @@ int win_gaussian(double *w, int n, int win_type, double alpha) switch(win_type & DSPL_WIN_SYM_MASK) { - case DSPL_WIN_SYMMETRIC: x = (double)(n-1)*0.5; break; - case DSPL_WIN_PERIODIC : x = (double)(n)*0.5; break; + case DSPL_WIN_SYMMETRIC: a = (double)(n-1)*0.5; break; + case DSPL_WIN_PERIODIC : a = (double)(n)*0.5; break; default: return ERROR_WIN_SYM; } - sigma = alpha / x; + sigma = 1.0 / (alpha * a); for(i = 0; i +#include +#include +#include + +#include "dspl.h" + +#define N 64 +#define K 1024 + + + +void win_norm(double* w, int n) +{ + int k; + double s = 0.0; + for(k = 0; k < n; k++) + s+=w[k]; + for(k = 0; k < n; k++) + w[k] /= s; +} + + + +void window_plot(double* w, int n, int k, int argc, char* argv[], + double ymin, double ymax, + char* png_fn, char* time_fn, char* freq_fn, char* title) +{ + void* hplot; /* GNUPLOT handles */ + char str[128] = {0}; + double* t = NULL; + double* W = NULL; + double* f = NULL; + double fs = (double)n; + fft_t pfft = {0}; + + t = (double*)malloc(n*sizeof(double)); + W = (double*)malloc(k*sizeof(double)); + f = (double*)malloc(k*sizeof(double)); + + linspace(0, n, n, DSPL_PERIODIC, t); + writetxt(t, w, n, time_fn); + win_norm(w, n); + fft_mag(w, k, &pfft, fs, DSPL_FLAG_LOGMAG | DSPL_FLAG_FFT_SHIFT, W, f); + writetxt(f, W, k, freq_fn); + + + /* plotting 3d surface by GNUPLOT */ + /* Create window 0 */ + gnuplot_create(argc, argv, 820, 360, png_fn, &hplot); + + gnuplot_cmd(hplot, "set grid"); + + memset(str,0,128); + sprintf(str, "set title '%s'", title); + gnuplot_cmd(hplot, str); + + gnuplot_cmd(hplot, "set multiplot layout 1,2 rowsfirst"); + gnuplot_cmd(hplot, "set xlabel 'n'"); + gnuplot_cmd(hplot, "set ylabel 'w(n)'"); + gnuplot_cmd(hplot, "unset key"); + + gnuplot_cmd(hplot, "set xrange[0:63]"); + + memset(str,0,128); + sprintf(str, "set yrange[%f:%f]", ymin, ymax); + gnuplot_cmd(hplot, str); + + memset(str,0,128); + sprintf(str, "plot '%s' w i lc 1,'%s' w p pt 7 ps 0.5 lc 1", + time_fn, time_fn); + gnuplot_cmd(hplot, str); + + gnuplot_cmd(hplot, "set xrange[-20:20]"); + gnuplot_cmd(hplot, "set yrange[-130:5]"); + memset(str,0,128); + + gnuplot_cmd(hplot, "set xlabel 'freq [DFT bins]'"); + gnuplot_cmd(hplot, "set ylabel 'W(freq), dB'"); + sprintf(str, "plot '%s' w l lc 2 ", freq_fn); + gnuplot_cmd(hplot, str); + gnuplot_close(hplot); + + if(t) + free(t); + if(W) + free(W); + if(f) + free(f); +} + + + +int main(int argc, char* argv[]) +{ + void* hdspl; /* DSPL handle */ + + hdspl = dspl_load(); /* Load DSPL function */ + double w[K] = {0}; + + /* Rectangular window */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_RECT, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_rect.png", + "dat/win_rect_time.txt", + "dat/win_rect_freq.txt", + "Rectangular window"); + + /* Bartlett window (triangular)*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BARTLETT, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_bartlett.png", + "dat/win_bartlett_time.txt", + "dat/win_bartlett_freq.txt", + "Bartlett window"); + + /* Flat top window */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_FLAT_TOP, 0.0); + window_plot(w, N, K, argc, argv, -1.0, 5.0, + "img/win_flattop.png", + "dat/win_flattop_time.txt", + "dat/win_flattop_freq.txt", + "Flat top window"); + + /* Bartlett - Hann window*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BARTLETT_HANN, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_bartletthann.png", + "dat/win_bartletthann_time.txt", + "dat/win_bartletthann_freq.txt", + "Bartlett-Hann window"); + + /* Blackman window*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BLACKMAN, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_blackman.png", + "dat/win_blackman_time.txt", + "dat/win_blackman_freq.txt", + "Blackman window"); + + /* Blackman - Harris window*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BLACKMAN_HARRIS, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_blackmanharris.png", + "dat/win_blackmanharris_time.txt", + "dat/win_blackmanharris_freq.txt", + "Blackman-Harris window"); + + + /* Blackman - Nuttall window*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BLACKMAN_NUTTALL, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_blackmannuttall.png", + "dat/win_blackmannuttall_time.txt", + "dat/win_blackmannuttall_freq.txt", + "Blackman-Nuttull window"); + + /* Dolph-Chebyshev window (sidelobes level is -50 dB)*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_CHEBY, 50.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_cheby50.png", + "dat/win_cheby50_time.txt", + "dat/win_cheby50_freq.txt", + "Dolph-Chebyshev window (Rs = 50dB)"); + + /* Dolph-Chebyshev window (sidelobes level is -80 dB)*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_CHEBY, 80.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_cheby80.png", + "dat/win_cheby80_time.txt", + "dat/win_cheby80_freq.txt", + "Dolph-Chebyshev window (Rs = 80dB)"); + + /* Dolph-Chebyshev window (sidelobes level is -120 dB)*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_CHEBY, 120.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_cheby120.png", + "dat/win_cheby120_time.txt", + "dat/win_cheby120_freq.txt", + "Dolph-Chebyshev window (Rs = 120dB)"); + + /* Gaussian window (sigma = 0.5)*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_GAUSSIAN, 0.5); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_gaussian0p5.png", + "dat/win_gaussian0p5_time.txt", + "dat/win_gaussian0p5_freq.txt", + "Gaussian window (sigma = 0.5)"); + + /* Gaussian window (sigma = 0.3)*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_GAUSSIAN, 0.3); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_gaussian0p3.png", + "dat/win_gaussian0p3_time.txt", + "dat/win_gaussian0p3_freq.txt", + "Gaussian window (sigma = 0.3)"); + + /* Hamming window*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_HAMMING, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_hamming.png", + "dat/win_hamming_time.txt", + "dat/win_hamming_freq.txt", + "Hamming window"); + + /* Hann window*/ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_HANN, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_hann.png", + "dat/win_hann_time.txt", + "dat/win_hann_freq.txt", + "Hann window"); + + /* Lanczos window */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_LANCZOS, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_lanczos.png", + "dat/win_lanczos_time.txt", + "dat/win_lanczos_freq.txt", + "Lanczos window"); + + /* Nuttall window */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_NUTTALL, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_nuttall.png", + "dat/win_nuttall_time.txt", + "dat/win_nuttall_freq.txt", + "Nuttall window"); + + /* Cosine window */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_COS, 0.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_cos.png", + "dat/win_cos_time.txt", + "dat/win_cos_freq.txt", + "Cosine window"); + + + /* Kaiser window pi * alpha = 4 */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_KAISER, 4.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_kaiser4p0.png", + "dat/win_kaiser4p0_time.txt", + "dat/win_kaiser4p0_freq.txt", + "Kaiser window (pi * alpha = 4)"); + + /* Kaiser window pi * alpha = 8 */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_KAISER, 8.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_kaiser8p0.png", + "dat/win_kaiser8p0_time.txt", + "dat/win_kaiser8p0_freq.txt", + "Kaiser window (pi * alpha = 8)"); + + /* Kaiser window pi * alpha = 12 */ + window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_KAISER, 12.0); + window_plot(w, N, K, argc, argv, 0.0, 1.1, + "img/win_kaiser12p0.png", + "dat/win_kaiser12p0_time.txt", + "dat/win_kaiser12p0_freq.txt", + "Kaiser window (pi * alpha = 12)"); + + + dspl_free(hdspl); /* free dspl handle */ + return 0; +} \ No newline at end of file