diff --git a/dspl/src/psd.c b/dspl/src/psd.c index f5200a9..60f1d84 100644 --- a/dspl/src/psd.c +++ b/dspl/src/psd.c @@ -33,7 +33,118 @@ #endif #ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup PSD_GROUP +\fn int psd_bartlett(double* x, int n, int nfft, + fft_t* pfft, double fs, + int flag, double* ppsd, double* pfrq) +\brief Непараметрическая оценка спектральной плотности мощности (СПМ) +вещественного сигнала методом Бартлетта. +Функция рассчитывает спектральную плотность мощности \f$ X(f) \f$ +выборки сигнала длительности \$n \$ отсчетов методом Бартлетта: +\f[ + X(f) = \frac{1}{N F_s } \sum_{p = 0}^{P-1}\left| \sum_{m = 0}^{n_{FFT}-1} + x(m+p \cdot n_{\text{FFT}}) \exp + \left( -j 2\pi f m \right) \right|^2, +\f] +где \f$ F_s \f$ -- частота +дискретизации (Гц), \f$P = n/n_{\text{FFT}}\f$ -- количество сегментов +смещений выборки сигналов размера \f$n_{FFT}\f$. + + +При использовании \f$n_{FFT} = n\f$ оценка Бартлетта переходит +в стандартную периодограмму. + +Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого +преобразования Фурье, для дискретной сетки частот от 0 Гц до \f$ F_s \f$ Гц +(по умолчанию), или от \f$-F_s /2 \f$ до \f$F_s /2 \f$, если установлен флаг +расчета двусторонней СПМ. + +\note Метод Бартлетта возвращает асимптотически несмещенную, +состоятельную оценку СПМ (уровень флуктуаций шумовой СПМ + уменьшается с ростом длины выборки `n` при фиксированной `nfft`). + +\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] nfft +Размер сегмента.\n +Размер выходного вектора СПМ, и соответствующего ей вектора частоты.\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_bartlett_test.c + +Программа производит расчет СПМ сигнала, состоящего из двух гармоник на +фоне белого гауссова шума. Расчет ведется по выборкe длины 8192 отсчета +при длине сегмента `nfft` 128, 1024 и +8192 отсчетов. + +Рассчитанные СПМ выводятся на графики: + +`n = 8192, nfft = 8192`: +\image html psd_bartlett_8192.png + +`n = 8192, nfft = 1024`: +\image html psd_bartlett_1024.png + +`n = 8192, nfft = 128`: +\image html psd_bartlett_1024.png + + +Можно видеть, что метод Бартлетта позволяет снизить +уровень флуктуация шума с увеличением количества сегментов. +Однако наблюдается эффект растекания спектра, который существенно ухудшает +динамический диапазон анализа. + +Для более качественной оценки СПМ смотри функцию \ref psd_welch + +\author Бахурин Сергей www.dsplib.org +***************************************************************************** */ #endif int DSPL_API psd_bartlett(double* x, int n, int nfft, fft_t* pfft, double fs, @@ -815,7 +926,7 @@ DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2) Пример периодограммных оценок СПМ для различной длины выборки сигнала: -\include psd_welch_cmplx_test.c +\include psd_welch_test.c Программа производит расчет СПМ сигнала, состоящего из двух комплексных гармоник на фоне белого гауссова шума. @@ -824,13 +935,13 @@ DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2) Рассчитанные СПМ выводятся на графики: `nfft = 8192, noverlap = 4096`: -\image html psd_welch_cmplx_8192.png +\image html psd_welch_8192.png `nfft = 1024, noverlap = 512`: -\image html psd_welch_cmplx_1024.png +\image html psd_welch_1024.png `nfft = 256, noverlap = 128`: -\image html psd_welch_cmplx_256.png +\image html psd_welch_256.png Можно видеть, что уменьшение `nfft` при фиксированной длительности сигнала позволяет уменьшить флуктуации шума и делает оценку состоятельной.