From 683f419847e222c5e5bfec1e27e79e551b9bcd64 Mon Sep 17 00:00:00 2001 From: Dsplib Date: Mon, 12 Oct 2020 22:40:56 +0300 Subject: [PATCH] added psd_welch example and doc Changes to be committed: modified: dspl/src/psd.c new file: examples/src/psd_welch_cmplx_test.c modified: examples/src/psd_welch_test.c --- dspl/src/psd.c | 609 ++++++++++++++++++++-------- examples/src/psd_welch_cmplx_test.c | 98 +++++ examples/src/psd_welch_test.c | 98 ++--- 3 files changed, 579 insertions(+), 226 deletions(-) create mode 100644 examples/src/psd_welch_cmplx_test.c diff --git a/dspl/src/psd.c b/dspl/src/psd.c index b80b93a..f5200a9 100644 --- a/dspl/src/psd.c +++ b/dspl/src/psd.c @@ -249,9 +249,9 @@ exit_label: /*! **************************************************************************** \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) + int win_type, double win_param, + fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) \brief Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом модифицированной периодограммы. @@ -709,38 +709,175 @@ exit_label: #endif #ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup PSD_GROUP +int psd_welch(double* x, int n, + int win_type, double win_param, + int nfft, int noverlap, fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) +\brief Непараметрическая оценка спектральной плотности мощности (СПМ) +вещественного сигнала методом Уэлча. + +Функция рассчитывает спектральную плотность мощности \f$ X(f) \f$ +выборки сигнала длительности \$n \$ отсчетов методом Уэлча: +\f[ + X(f) = \frac{1}{U P F_s } \sum_{p = 0}^{P-1}\left| \sum_{m = 0}^{n_{FFT}-1} + w(m) x(m+p \cdot n_{\text{overlap}}) \exp + \left( -j 2\pi f m \right) \right|^2, +\f] +где \f$ w(m) \f$ -- отсчёты оконной функции, \f$ F_s \f$ -- частота +дискретизации (Гц), \f$P = n/n_{\text{overlap}}\f$ -- количество сегментов +смещений выборки сигналов размера \f$n_{FFT}\f$, + + \f$ U \f$ нормировочный коэффициент равный +\f[ + U = \sum_{m = 0}^{n-1} w^2(m), +\f] + +Процедура разбиения исходной последовательности длительности `n` отсчетов +на сегменты длины \f$n_{FFT}\f$ отсчетов, перекрывающихся с интервалом +\f$n_{\text{overlap}}\f$ отсчетов, показан на следующем рисунке + +\image html welch_overlap.png + +Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого +преобразования Фурье, для дискретной сетки частот от 0 Гц до \f$ F_s \f$ Гц +(по умолчанию), или от \f$-F_s /2 \f$ до \f$F_s /2 \f$, если установлен флаг +расчета двусторонней периодограммы. + +\note Периодограмма Уэлча возвращает смещенную, но состоятельную оценку СПМ. + +\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] nfft +Размер перекрывающегося сегмента.\n +Размер выходного вектора СПМ, и соответсвующего ей вектора частоты.\n\n + +\param[in] noverlap +Размер сдвига сегментов относительно друг друга (отсчетов).\n +`noverlap = nfft` задает оценку без перекрытия сегментов. \n +Обычно используют сдвиг равный половине размера сегмента `noverlap = nfft/2`.\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 +Размер вектора `[nfft x 1]`. \n +Память должна быть выделена. \n\n + +\param[in, out] pfrq +Указатель на вектор частоты, соответствующей +значениям рассчитанного вектора СПМ. \n +Размер вектора `[nfft x 1]`. \n +Указатель может быть `NULL`,в этом случае вектор частоты не +рассчитывается и не возвращается. \n\n + +\return +`RES_OK` если расчет произведен успешно. \n +В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n \n + +Пример периодограммных оценок СПМ для различной длины выборки сигнала: + +\include psd_welch_cmplx_test.c + +Программа производит расчет СПМ сигнала, состоящего из двух комплексных +гармоник на фоне белого гауссова шума. +Расчет ведется по выборке сигнала длины 8192 отсчета. + +Рассчитанные СПМ выводятся на графики: + +`nfft = 8192, noverlap = 4096`: +\image html psd_welch_cmplx_8192.png + +`nfft = 1024, noverlap = 512`: +\image html psd_welch_cmplx_1024.png + +`nfft = 256, noverlap = 128`: +\image html psd_welch_cmplx_256.png + +Можно видеть, что уменьшение `nfft` при фиксированной длительности сигнала +позволяет уменьшить флуктуации шума и делает оценку состоятельной. + +\author Бахурин Сергей www.dsplib.org +***************************************************************************** */ #endif int DSPL_API psd_welch(double* x, int n, int win_type, double win_param, - int npsd, int noverlap, fft_t* pfft, double fs, + int nfft, int noverlap, fft_t* pfft, double fs, int flag, double* ppsd, double* pfrq) { - double *win = NULL; - double wg; - complex_t *tmp = NULL; - double *s = NULL; - int err, k, pos, cnt; + int err, pos, cnt, k; + double *pdgr = NULL; + double *tmp = NULL; + double *w = NULL; fft_t *ptr_fft = NULL; + double wn; + + pos = cnt = 0; + + pdgr = (double*)malloc(nfft * sizeof(double)); + if(!pdgr) + return ERROR_MALLOC; + + tmp = (double*) malloc(nfft * sizeof(double)); + if(!tmp) + return ERROR_MALLOC; - if(!x || !ppsd) - return ERROR_PTR; - - if(n<1 || npsd < 1) - return ERROR_SIZE; - - if(noverlap < 1 || noverlap > npsd) - return ERROR_OVERLAP; - - if(fs < 0.0) - return ERROR_FS; - - win = (double*)malloc(npsd * sizeof(double)); - if(!win) + /* window malloc */ + w = (double*)malloc(nfft*sizeof(double)); + if(!w) { err = ERROR_MALLOC; goto exit_label; } + + /* create window */ + err = window(w, nfft, win_type, win_param); + if(err != RES_OK) + goto exit_label; + + /* window normalization wn = sum(w.^2) */ + wn = 0.0; + for(k = 0; k < nfft; k++) + wn += w[k]*w[k]; + if(!pfft) { @@ -749,142 +886,261 @@ int DSPL_API psd_welch(double* x, int n, } else ptr_fft = pfft; - - - err = window(win, npsd, win_type, win_param); - if(err != RES_OK) - goto exit_label; - - wg = 0.0; - for(k = 0; k < npsd; k++) - wg += win[k] * win[k]; - wg = 1.0 / wg; - tmp = (complex_t*)malloc(npsd*sizeof(complex_t)); - if(!tmp) + + + memset(ppsd, 0, nfft * sizeof(double)); + while(pos + nfft <= n) { - err = ERROR_MALLOC; - goto exit_label; - } - - s = (double*)malloc(npsd*sizeof(double)); - if(!s) - { - err = ERROR_MALLOC; - goto exit_label; - } - - pos = 0; - cnt = 0; - memset(ppsd, 0, npsd*sizeof(double)); - while(pos+npsd <= n) - { - for(k = 0; k < npsd; k++) - s[k] = x[k+pos]*win[k]; - - err = fft(s, npsd, ptr_fft, tmp); + for(k = 0; k < nfft; k++) + tmp[k] = x[pos+k] * w[k]; + err = fft_mag(tmp, nfft, ptr_fft, fs, + flag & DSPL_FLAG_FFT_SHIFT, pdgr, NULL); if(err != RES_OK) goto exit_label; - - for(k = 0; k < npsd; k++) - ppsd[k] += wg * ABSSQR(tmp[k]); - + for(k = 0; k < nfft; k++) + ppsd[k] += pdgr[k]; pos += noverlap; cnt++; - - } - - for(k = 0; k < npsd; k++) - ppsd[k] /= (double)cnt * fs; - - if(flag & DSPL_FLAG_LOGMAG) - { - for(k = 0; k < npsd; k++) - ppsd[k] = 10.0 * log10(ppsd[k] + DBL_EPSILON); } - - - if(flag & DSPL_FLAG_PSD_TWOSIDED) + if(pos < n) { - err = fft_shift(ppsd, npsd, ppsd); + + memset(tmp ,0, nfft * sizeof(double)); + for(k = 0; k < n - pos; k++) + tmp[k] = x[pos+k] * w[k]; + + + err = fft_mag(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]; + + cnt++; } + /* fill frequency */ if(pfrq) { - if(flag & DSPL_FLAG_PSD_TWOSIDED) - { - err = linspace(-0.5*fs, fs*0.5, npsd, DSPL_PERIODIC, pfrq); - if(err != RES_OK) - goto exit_label; - } + 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, npsd, DSPL_PERIODIC, pfrq); - if(err != RES_OK) - goto exit_label; - } + err = linspace(0, fs, nfft, DSPL_PERIODIC, pfrq); + } + + /* scale magnitude */ + if(flag & DSPL_FLAG_LOGMAG) + { + printf("wn = %f\n", wn); + for(k = 0; k < nfft; k++) + ppsd[k] = 10.0 * log10(ppsd[k] / (fs * wn * (double)cnt)); + } + else + { + for(k = 0; k < nfft; k++) + ppsd[k] /= fs * wn * (double)cnt; } - err = RES_OK; + exit_label: - if(win) - free(win); + if(pdgr) + free(pdgr); if(tmp) free(tmp); - if(s) - free(s); + if(w) + free(w); 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 +int psd_welch_cmplx(complex_t* x, int n, + int win_type, double win_param, + int nfft, int noverlap, fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) +\brief Непараметрическая оценка спектральной плотности мощности (СПМ) +комплексного сигнала методом Уэлча. +Функция рассчитывает спектральную плотность мощности \f$ X(f) \f$ +выборки сигнала длительности \$n \$ отсчетов методом Уэлча: +\f[ + X(f) = \frac{1}{U P F_s } \sum_{p = 0}^{P-1}\left| \sum_{m = 0}^{n_{FFT}-1} + w(m) x(m+p \cdot n_{\text{overlap}}) \exp + \left( -j 2\pi f m \right) \right|^2, +\f] +где \f$ w(m) \f$ -- отсчёты оконной функции, \f$ F_s \f$ -- частота +дискретизации (Гц), \f$P = n/n_{\text{overlap}}\f$ -- количество сегментов +смещений выборки сигналов размера \f$n_{FFT}\f$, + + \f$ U \f$ нормировочный коэффициент равный +\f[ + U = \sum_{m = 0}^{n-1} w^2(m), +\f] + +Процедура разбиения исходной последовательности длительности `n` отсчетов +на сегменты длины \f$n_{FFT}\f$ отсчетов, перекрывающихся с интервалом +\f$n_{\text{overlap}}\f$ отсчетов, показан на следующем рисунке + +\image html welch_overlap.png + +Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого +преобразования Фурье, для дискретной сетки частот от 0 Гц до \f$ F_s \f$ Гц +(по умолчанию), или от \f$-F_s /2 \f$ до \f$F_s /2 \f$, если установлен флаг +расчета двусторонней периодограммы. + +\note Периодограмма Уэлча возвращает смещенную, но состоятельную оценку СПМ. + +\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] nfft +Размер перекрывающегося сегмента.\n +Размер выходного вектора СПМ, и соответсвующего ей вектора частоты.\n\n + +\param[in] noverlap +Размер сдвига сегментов относительно друг друга (отсчетов).\n +`noverlap = nfft` задает оценку без перекрытия сегментов. \n +Обычно используют сдвиг равный половине размера сегмента `noverlap = nfft/2`.\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 +Размер вектора `[nfft x 1]`. \n +Память должна быть выделена. \n\n + +\param[in, out] pfrq +Указатель на вектор частоты, соответствующей +значениям рассчитанного вектора СПМ. \n +Размер вектора `[nfft x 1]`. \n +Указатель может быть `NULL`,в этом случае вектор частоты не +рассчитывается и не возвращается. \n\n + +\return +`RES_OK` если расчет произведен успешно. \n +В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n \n + +Пример периодограммных оценок СПМ для различной длины выборки сигнала: + +\include psd_welch_cmplx_test.c + +Программа производит расчет СПМ сигнала, состоящего из двух комплексных +гармоник на фоне белого гауссова шума. +Расчет ведется по выборке сигнала длины 8192 отсчета. + +Рассчитанные СПМ выводятся на графики: + +`nfft = 8192, noverlap = 4096`: +\image html psd_welch_cmplx_8192.png + +`nfft = 1024, noverlap = 512`: +\image html psd_welch_cmplx_1024.png + +`nfft = 256, noverlap = 128`: +\image html psd_welch_cmplx_256.png + +Можно видеть, что уменьшение `nfft` при фиксированной длительности сигнала +позволяет уменьшить флуктуации шума и делает оценку состоятельной. + +\author Бахурин Сергей www.dsplib.org +***************************************************************************** */ #endif int DSPL_API psd_welch_cmplx(complex_t* x, int n, - int win_type, double win_param, - int npsd, int noverlap, fft_t* pfft, double fs, - int flag, double* ppsd, double* pfrq) + int win_type, double win_param, + int nfft, int noverlap, fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) { - double *win = NULL; - complex_t *tmp = NULL; - complex_t *s = NULL; - fft_t *ptr_fft = NULL; + int err, pos, cnt, k; + double *pdgr = NULL; + complex_t *tmp = NULL; + double *w = NULL; + fft_t *ptr_fft = NULL; + double wn; - double wg; - int err, k, pos, cnt; + pos = cnt = 0; - if(!x || !ppsd) - return ERROR_PTR; + pdgr = (double*)malloc(nfft * sizeof(double)); + if(!pdgr) + return ERROR_MALLOC; + + tmp = (complex_t*) malloc(nfft * sizeof(complex_t)); + if(!tmp) + return ERROR_MALLOC; - if(n<1 || npsd < 1) - return ERROR_SIZE; - if(noverlap < 1 || noverlap > npsd) - return ERROR_OVERLAP; - - if(fs < 0.0) - return ERROR_FS; - - win = (double*)malloc(npsd * sizeof(double)); - if(!win) + /* window malloc */ + w = (double*)malloc(nfft*sizeof(double)); + if(!w) { err = ERROR_MALLOC; goto exit_label; } + + /* create window */ + err = window(w, nfft, win_type, win_param); + if(err != RES_OK) + goto exit_label; + + /* window normalization wn = sum(w.^2) */ + wn = 0.0; + for(k = 0; k < nfft; k++) + wn += w[k]*w[k]; + if(!pfft) { @@ -893,93 +1149,88 @@ int DSPL_API psd_welch_cmplx(complex_t* x, int n, } else ptr_fft = pfft; - - err = window(win, npsd, win_type, win_param); - if(err != RES_OK) - goto exit_label; - - wg = 0.0; - for(k = 0; k < npsd; k++) - wg += win[k] * win[k]; - wg = 1.0 / wg; - - tmp = (complex_t*)malloc(npsd*sizeof(complex_t)); - if(!tmp) + + + + + memset(ppsd, 0, nfft * sizeof(double)); + while(pos + nfft <= n) { - err = ERROR_MALLOC; - goto exit_label; - } - - s = (complex_t*)malloc(npsd*sizeof(complex_t)); - if(!s) - { - err = ERROR_MALLOC; - goto exit_label; - } - - pos = 0; - cnt = 0; - memset(ppsd, 0, npsd*sizeof(double)); - while(pos+npsd <= n) - { - for(k = 0; k < npsd; k++) + for(k = 0; k < nfft; k++) { - RE(s[k]) = RE(x[k+pos])*win[k]; - IM(s[k]) = IM(x[k+pos])*win[k]; + RE(tmp[k]) = RE(x[pos+k]) * w[k]; + IM(tmp[k]) = IM(x[pos+k]) * w[k]; } - err = fft_cmplx(s, npsd, ptr_fft, tmp); + 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 < npsd; k++) - ppsd[k] += wg * ABSSQR(tmp[k]); - + for(k = 0; k < nfft; k++) + ppsd[k] += pdgr[k]; pos += noverlap; cnt++; - } - - for(k = 0; k < npsd; k++) - ppsd[k] /= (double)cnt * fs; - - if(flag & DSPL_FLAG_LOGMAG) + + if(pos < n) { - for(k = 0; k < npsd; k++) - ppsd[k] = 10.0 * log10(ppsd[k] + DBL_EPSILON); - } - if(flag & DSPL_FLAG_PSD_TWOSIDED) - { - err = fft_shift(ppsd, npsd, ppsd); + memset(tmp ,0, nfft * sizeof(complex_t)); + for(k = 0; k < n - pos; k++) + { + RE(tmp[k]) = RE(x[pos+k]) * w[k]; + IM(tmp[k]) = IM(x[pos+k]) * w[k]; + } + + 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]; + + cnt++; } + /* fill frequency */ if(pfrq) { - if(flag & DSPL_FLAG_PSD_TWOSIDED) - { - err = linspace(-0.5*fs, fs*0.5, npsd, DSPL_PERIODIC, pfrq); - if(err != RES_OK) - goto exit_label; - } + 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, npsd, DSPL_PERIODIC, pfrq); - if(err != RES_OK) - goto exit_label; - } + err = linspace(0, fs, nfft, DSPL_PERIODIC, pfrq); } - err = RES_OK; + + /* scale magnitude */ + if(flag & DSPL_FLAG_LOGMAG) + { + for(k = 0; k < nfft; k++) + ppsd[k] = 10.0 * log10(ppsd[k] / (fs * wn * (double)cnt)); + } + else + { + for(k = 0; k < nfft; k++) + ppsd[k] /= fs * wn * (double)cnt; + } + exit_label: - if(win) - free(win); + if(pdgr) + free(pdgr); if(tmp) free(tmp); - if(s) - free(s); + if(w) + free(w); if(ptr_fft && (ptr_fft != pfft)) + { + fft_free(ptr_fft); free(ptr_fft); + } return err; } + diff --git a/examples/src/psd_welch_cmplx_test.c b/examples/src/psd_welch_cmplx_test.c new file mode 100644 index 0000000..4213fbb --- /dev/null +++ b/examples/src/psd_welch_cmplx_test.c @@ -0,0 +1,98 @@ +#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.20 * (double)k); + + IM(x[k]) += sin(M_2PI * 0.26 * (double)k) + + 0.1 * sin(M_2PI * 0.20 * (double)k); + } + + /* Twosided Welch PSD logscale magnitude + n = 8192, nfft = 8192, noverlap = 4096 */ + err = psd_welch_cmplx(x, N, DSPL_WIN_BLACKMAN, 0, 8192, 4096, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, 8192, "dat/psd_welch_cmplx_8192.txt"); + psd_plot(argc, argv, "img/psd_welch_cmplx_8192.png", + "dat/psd_welch_cmplx_8192.txt"); + + + /* Twosided Welch PSD logscale magnitude + n = 8192, nfft = 1024, noverlap = 512 */ + err = psd_welch_cmplx(x, N, DSPL_WIN_BLACKMAN, 0, 1024, 512, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, 1024, "dat/psd_welch_cmplx_1024.txt"); + psd_plot(argc, argv, "img/psd_welch_cmplx_1024.png", + "dat/psd_welch_cmplx_1024.txt"); + + + /* Twosided Welch PSD logscale magnitude + n = 8192, nfft = 256, noverlap = 128 */ + err = psd_welch_cmplx(x, N, DSPL_WIN_BLACKMAN, 0, 256, 128, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, 256, "dat/psd_welch_cmplx_256.txt"); + psd_plot(argc, argv, "img/psd_welch_cmplx_256.png", + "dat/psd_welch_cmplx_256.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_welch_test.c b/examples/src/psd_welch_test.c index 5137b6f..9899bb9 100644 --- a/examples/src/psd_welch_test.c +++ b/examples/src/psd_welch_test.c @@ -4,80 +4,84 @@ #include "dspl.h" -#define N 8192 -#define NFFT 512 +#define N 8192 +#define FS 1.0 -#define FS 0.01 +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 */ - void* hplot; /* GNUPLOT handle */ - hdspl = dspl_load(); // Load DSPL function + void* hdspl; /* DSPL handle */ + + hdspl = dspl_load(); /* Load DSPL function */ random_t rnd = {0}; /* random structure */ - double *x = NULL; + double *x = NULL; double *psd = NULL; double *frq = NULL; - + int k, err; - x = (double*)malloc(N*sizeof(double)); + 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 */ + /* 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.1*k) + cos(2*pi*0.2*k) + noise */ + /* 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); + x[k] += cos(M_2PI * 0.26 * (double)k) + 0.1*cos(M_2PI*0.2* (double)k); - /* Twosided PSD with logscale magnitude */ - err = psd_welch(x, N, DSPL_WIN_BLACKMAN , 0, NFFT, NFFT, NULL, FS, - DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); - - printf("error: 0x%.8x\n", err); - - // Save PSD to the "dat/psd_welch.txt" file - writetxt(frq, psd, NFFT, "dat/psd_welch.txt"); - - - /* Twosided PSD with logscale magnitude */ - err = psd_periodogram(x, N, DSPL_WIN_BLACKMAN, 0, NULL, FS, - DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); - printf("error: 0x%.8x\n", err); - writetxt(frq, psd, N, "dat/psd_periodogram.txt"); - - - /* Twosided PSD with logscale magnitude */ - err = psd_bartlett(x, N, NFFT, NULL, FS, + /* Twosided PSD logscale magnitude n = 8192, nfft = 8192. + This case is periodogram */ + + err = psd_welch(x, N, DSPL_WIN_BLACKMAN, 0, 8192, 4096, NULL, FS, DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); - printf("error: 0x%.8x\n", err); + writetxt(frq, psd, 8192, "dat/psd_welch_8192.txt"); + psd_plot(argc, argv, "img/psd_welch_8192.png", + "dat/psd_welch_8192.txt"); - // Save PSD to the "dat/psd_welch.txt" file - writetxt(frq, psd, NFFT, "dat/psd_barlett.txt"); - /* plotting by GNUPLOT */ - gnuplot_create(argc, argv, 560, 420, "img/psd_welch.png", &hplot); - gnuplot_cmd(hplot, "unset key"); - gnuplot_cmd(hplot, "set grid"); - gnuplot_cmd(hplot, "set xlabel 'frequency'"); - gnuplot_cmd(hplot, "set ylabel 'PSD, dB'"); + err = psd_welch(x, N, DSPL_WIN_BLACKMAN, 0, 1024, 512, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, 1024, "dat/psd_welch_1024.txt"); + psd_plot(argc, argv, "img/psd_welch_1024.png", + "dat/psd_welch_1024.txt"); - gnuplot_cmd(hplot, "plot 'dat/psd_periodogram.txt' w l lt 3,\\"); - gnuplot_cmd(hplot, " 'dat/psd_barlett.txt' w l lt 1,\\"); - gnuplot_cmd(hplot, " 'dat/psd_welch.txt' w l lt -1"); + err = psd_welch(x, N, DSPL_WIN_BLACKMAN, 0, 256, 128, NULL, FS, + DSPL_FLAG_LOGMAG | DSPL_FLAG_PSD_TWOSIDED, psd, frq); + writetxt(frq, psd, 256, "dat/psd_welch_256.txt"); + psd_plot(argc, argv, "img/psd_welch_256.png", + "dat/psd_welch_256.txt"); - gnuplot_close(hplot); - - dspl_free(hdspl); // free dspl handle + dspl_free(hdspl); /* free dspl handle */ if(x) free(x); - + if(psd) + free(psd); + if(frq) + free(frq); return 0; }