Функции
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)
	Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом Уэлча. Подробнее...

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)
	Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом Уэлча. Подробнее...

int	psd_bartlett (double x, int n, int nfft, fft_t pfft, double fs, int flag, double ppsd, double pfrq)
	Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом Бартлетта. Подробнее...

int	psd_bartlett_cmplx (complex_t x, int n, int nfft, fft_t pfft, double fs, int flag, double ppsd, double pfrq)
	Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом Бартлетта. Подробнее...

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	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)
	Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом модифицированной периодограммы. Подробнее...

Подробное описание

Функции

◆ psd_bartlett()

int psd_bartlett	(	double *	x,
		int	n,
		int	nfft,
		fft_t *	pfft,
		double	fs,
		int	flag,
		double *	ppsd,
		double *	pfrq
	)

Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом Бартлетта.

Функция рассчитывает спектральную плотность мощности $ X(f) $ выборки сигнала длительности $n $ отсчетов методом Бартлетта:

\[ 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_s $ – частота дискретизации (Гц), $P = n/n_{\text{FFT}}$ – количество сегментов смещений выборки сигналов размера $n_{FFT}$.

При использовании $n_{FFT} = n$ оценка Бартлетта переходит в стандартную периодограмму.

Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого преобразования Фурье, для дискретной сетки частот от 0 Гц до $ F_s $ Гц (по умолчанию), или от $-F_s /2 $ до $F_s /2 $, если установлен флаг расчета двусторонней СПМ.

Заметки: Метод Бартлетта возвращает асимптотически несмещенную, состоятельную оценку СПМ (уровень флуктуаций шумовой СПМ уменьшается с ростом длины выборки n при фиксированной nfft).

Аргументы

[in]	x	Указатель на входной вектор вещественного сигнала $x(m)$, $ m = 0 \ldots n-1 $. Размер вектора `[n x 1]`.
[in]	n	Размер вектора входного сигнала. Также размер выходного вектора СПМ и вектора частоты также равны `n`.
[in]	nfft	Размер сегмента. Размер выходного вектора СПМ, и соответствующего ей вектора частоты.
[in]	pfft	Указатель на структуру fft_t. Указатель может быть `NULL`. В этом случае объект структуры будет создан внутри функции и удален перед завершением. Если предполагается многократный вызов функции, то рекомендуется создать объект fft_t и передавать в функцию, чтобы не создавать его каждый раз.
[in]	fs	частота дискретизации выборки исходного сигнала (Гц).
[in]	flag	Комбинация битовых флагов, задающих режим расчета: DSPL_FLAG_LOGMAG - СПМ считать в логарифмическом масштабе в единицах дБ/Гц DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2)
[in,out]	ppsd	Указатель на вектор СПМ рассчитанных по входному сигналу $x$. Размер вектора `[n x 1]`. Память должна быть выделена.
[in,out]	pfrq	Указатель на вектор частоты, соответствующей значениям рассчитанного вектора СПМ. Размер вектора `[n x 1]`. Указатель может быть `NULL`,в этом случае вектор частоты не рассчитывается и не возвращается.

Возвращает: RES_OK если расчет произведен успешно.
В противном случае код ошибки.

Пример оценок СПМ методом Бартлетта:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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;
}
 
 

Программа производит расчет СПМ сигнала, состоящего из двух гармоник на фоне белого гауссова шума. Расчет ведется по выборкe длины 8192 отсчета при длине сегмента nfft 128, 1024 и 8192 отсчетов.

Рассчитанные СПМ выводятся на графики:

n = 8192, nfft = 8192:

n = 8192, nfft = 1024:

n = 8192, nfft = 128:

Можно видеть, что метод Бартлетта позволяет снизить уровень флуктуация шума с увеличением количества сегментов. Однако наблюдается эффект растекания спектра, который существенно ухудшает динамический диапазон анализа.

Для более качественной оценки СПМ смотри функцию psd_welch

Автор: Бахурин Сергей www.dsplib.org

См. определение в файле psd_bartlett.c строка 146

◆ psd_bartlett_cmplx()

int psd_bartlett_cmplx	(	complex_t *	x,
		int	n,
		int	nfft,
		fft_t *	pfft,
		double	fs,
		int	flag,
		double *	ppsd,
		double *	pfrq
	)

Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом Бартлетта.

Функция рассчитывает спектральную плотность мощности $ X(f) $ выборки сигнала длительности $n $ отсчетов методом Бартлетта:

\[ 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_s $ – частота дискретизации (Гц), $P = n/n_{\text{FFT}}$ – количество сегментов смещений выборки сигналов размера $n_{FFT}$.

При использовании $n_{FFT} = n$ оценка Бартлетта переходит в стандартную периодограмму.

Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого преобразования Фурье, для дискретной сетки частот от 0 Гц до $ F_s $ Гц (по умолчанию), или от $-F_s /2 $ до $F_s /2 $, если установлен флаг расчета двусторонней СПМ.

Заметки: Метод Бартлетта возвращает асимптотически несмещенную, состоятельную оценку СПМ (уровень флуктуаций шумовой СПМ уменьшается с ростом длины выборки n при фиксированной nfft).

Аргументы

[in]	x	Указатель на входной вектор комплексного сигнала $x(m)$, $ m = 0 \ldots n-1 $. Размер вектора `[n x 1]`.
[in]	n	Размер вектора входного сигнала. Также размер выходного вектора СПМ и вектора частоты также равны `n`.
[in]	nfft	Размер сегмента. Размер выходного вектора СПМ, и соответствующего ей вектора частоты.
[in]	pfft	Указатель на структуру fft_t. Указатель может быть `NULL`. В этом случае объект структуры будет создан внутри функции и удален перед завершением. Если предполагается многократный вызов функции, то рекомендуется создать объект fft_t и передавать в функцию, чтобы не создавать его каждый раз.
[in]	fs	частота дискретизации выборки исходного сигнала (Гц).
[in]	flag	Комбинация битовых флагов, задающих режим расчета: DSPL_FLAG_LOGMAG - СПМ считать в логарифмическом масштабе в единицах дБ/Гц DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2)
[in,out]	ppsd	Указатель на вектор СПМ рассчитанных по входному сигналу $x$. Размер вектора `[n x 1]`. Память должна быть выделена.
[in,out]	pfrq	Указатель на вектор частоты, соответствующей значениям рассчитанного вектора СПМ. Размер вектора `[n x 1]`. Указатель может быть `NULL`,в этом случае вектор частоты не рассчитывается и не возвращается.

Возвращает: RES_OK если расчет произведен успешно.
В противном случае код ошибки.

Пример оценок СПМ методом Бартлетта:

Программа производит расчет СПМ сигнала, состоящего из двух комплексных экспонент на фоне белого гауссова шума. Расчет ведется по выборкe длины 8192 отсчета при длине сегмента nfft 128, 1024 и 8192 отсчетов.

Рассчитанные СПМ выводятся на графики:

n = 8192, nfft = 8192:

n = 8192, nfft = 1024:

n = 8192, nfft = 128:

Можно видеть, что метод Бартлетта позволяет снизить уровень флуктуация шума с увеличением количества сегментов. Однако наблюдается эффект растекания спектра, который существенно ухудшает динамический диапазон анализа.

Для более качественной оценки СПМ смотри функцию psd_welch_cmplx

Автор: Бахурин Сергей www.dsplib.org

См. определение в файле psd_bartlett_cmplx.c строка 145

◆ psd_periodogram()

int psd_periodogram	(	double *	x,
		int	n,
		int	win_type,
		double	win_param,
		fft_t *	pfft,
		double	fs,
		int	flag,
		double *	ppsd,
		double *	pfrq
	)

Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом модифицированной периодограммы.

Функция рассчитывает спектральную плотность мощности $ X(f) $ выборки сигнала длительности $n $ отсчетов методом модифицированной периодограммы:

\[ 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, \]

где $ w(m) $ – отсчёты оконной функции, $ F_s $ – частота дискретизации (Гц), $ U $ нормировочный коэффициент равный

\[ U = \sum_{m = 0}^{n-1} w^2(m) \]

При использовании прямоугольного окна модифицированная периодограмма переходит в стандартную периодограмму.

Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого преобразования Фурье, для дискретной сетки частот от 0 Гц до $ F_s $ Гц (по умолчанию), или от $-F_s /2 $ до $F_s /2 $, если установлен флаг расчета двусторонней периодограммы.

Заметки: Периодограмма возвращает асимптотически несмещенную, но несостоятельную оценку СПМ (уровень флуктуаций шумовой составляющей СПМ не уменьшается с ростом длины выборки n).

Аргументы

[in]	x	Указатель на входной вектор вещественного сигнала $x(m)$, $ m = 0 \ldots n-1 $. Размер вектора `[n x 1]`.
[in]	n	Размер вектора входного сигнала. Также размер выходного вектора СПМ и вектора частоты также равны `n`.
[in]	win_type	Тип оконной функции, применяемой для модифицированной периодограммы. Подробнее смотри описание функции window.
[in]	win_param	Параметр оконной функции. Данный параметр применяется только для параметрических типов окон (смотри описание функции window). Для непараметрических функций игнорируется.
[in]	pfft	Указатель на структуру fft_t. Указатель может быть `NULL`. В этом случае объект структуры будет создан внутри функции и удален перед завершением. Если предполагается многократный вызов функции, то рекомендуется создать объект fft_t и передавать в функцию, чтобы не создавать его каждый раз.
[in]	fs	частота дискретизации выборки исходного сигнала (Гц).
[in]	flag	Комбинация битовых флагов, задающих режим расчета: DSPL_FLAG_LOGMAG - СПМ считать в логарифмическом масштабе в единицах дБ/Гц DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2)
[in,out]	ppsd	Указатель на вектор СПМ рассчитанных по входному сигналу $x$. Размер вектора `[n x 1]`. Память должна быть выделена.
[in,out]	pfrq	Указатель на вектор частоты, соответствующей значениям рассчитанного вектора СПМ. Размер вектора `[n x 1]`. Указатель может быть `NULL`,в этом случае вектор частоты не рассчитывается и не возвращается.

Возвращает: RES_OK если расчет произведен успешно.
В противном случае код ошибки.

Пример периодограммных оценок СПМ для различной длины выборки сигнала:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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;
}
 
 

Программа производит расчет СПМ сигнала, состоящего из двух гармоник на фоне белого гауссова шума. Расчет ведется по выборкам длины 128, 1024 и 8192 отсчетов.

В результате периодограммы (стандартная с прямоугольным окном и модифицированная с окном Блэкмана) выводятся на графики:

n = 8192 точек (черная – классическая периодограмма с прямоугольным окном, зеленая – модифицированная с окном Блэкмана):

n = 1024 точек (черная – классическая периодограмма с прямоугольным окном, зеленая – модифицированная с окном Блэкмана):

n = 128 точек (черная – классическая периодограмма с прямоугольным окном, зеленая – модифицированная с окном Блэкмана):

Можно видеть, что модифицированная периодограмма позволяет снизить растекание СПМ, однако уровень флуктуация шума не уменьшается с увеличением размера выборки от 128 до 8192 отсчетов (оценка несостоятельная).

Автор: Бахурин Сергей www.dsplib.org

См. определение в файле psd_periodogram.c строка 152

◆ psd_periodogram_cmplx()

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
	)

Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом модифицированной периодограммы.

Функция рассчитывает спектральную плотность мощности $ X(f) $ выборки сигнала длительности $n $ отсчетов методом модифицированной периодограммы:

\[ 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, \]

где $ w(m) $ – отсчёты оконной функции, $ F_s $ – частота дискретизации (Гц), $ U $ нормировочный коэффициент равный

\[ U = \sum_{m = 0}^{n-1} w^2(m) \]

При использовании прямоугольного окна модифицированная периодограмма переходит в стандартную периодограмму.

Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого преобразования Фурье, для дискретной сетки частот от 0 Гц до $ F_s $ Гц (по умолчанию), или от $-F_s /2 $ до $F_s /2 $, если установлен флаг расчета двусторонней периодограммы.

Заметки: Периодограмма возвращает асимптотически несмещенную, но несостоятельную оценку СПМ (уровень флуктуаций шумовой составляющей СПМ не уменьшается с ростом длины выборки n).

Аргументы

[in]	x	Указатель на входной вектор вещественного сигнала $x(m)$, $ m = 0 \ldots n-1 $. Размер вектора `[n x 1]`.
[in]	n	Размер вектора входного сигнала. Также размер выходного вектора СПМ и вектора частоты также равны `n`.
[in]	win_type	Тип оконной функции, применяемой для модифицированной периодограммы. Подробнее смотри описание функции window.
[in]	win_param	Параметр оконной функциии. Данный параметр применяется только для парамтрических типов окон (смотри описание функции window). Для непараметрических функций игнорируется.
[in]	pfft	Указатель на структуру fft_t. Указатель может быть `NULL`. В этом случае объект структуры будет создан внутри функции и удален перед завершением. Если предполагается многократный вызов функции, то рекомендуется создать объект fft_t и передавать в функцию, чтобы не создавать его каждый раз.
[in]	fs	частота дискретизации выборки исходного сигнала (Гц).
[in]	flag	Комбинация битовых флагов, задающих режим расчета: DSPL_FLAG_LOGMAG - СПМ считать в логарифмическом масштабе в единицах дБ/Гц DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2)
[in,out]	ppsd	Указатель на вектор СПМ рассчитанных по входному сигналу $x$. Размер вектора `[n x 1]`. Память должна быть выделена.
[in,out]	pfrq	Указатель на вектор частоты, соответствующей значениям рассчитанного вектора СПМ. Размер вектора `[n x 1]`. Указатель может быть `NULL`,в этом случае вектор частоты не рассчитывается и не возвращается.

Возвращает: RES_OK если расчет произведен успешно.
В противном случае код ошибки.

Пример периодограммных оценок СПМ для различной длины выборки сигнала:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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;
}
 
 

Программа производит расчет СПМ сигнала, состоящего из двух гармоник на фоне белого гауссова шума. Расчет ведется по выборкам длины 128, 1024 и 8192 отсчетов.

В результате периодограммы (стандартная с прямоугольным окном и модифицированная с окном Блэкмана) выводятся на графики:

n = 8192 точек (черная – классическая периодограмма с прямоугольным окном, зеленая – модифицированная с окном Блэкмана):

n = 1024 точек (черная – классическая периодограмма с прямоугольным окном, зеленая – модифицированная с окном Блэкмана):

n = 128 точек (черная – классическая периодограмма с прямоугольным окном, зеленая – модифицированная с окном Блэкмана):

Можно видеть, что модифицированная периодограмма позволяет снизить растекание СПМ, однако уровень флуктуация шума не уменьшается с увеличением размера выборки от 128 до 8192 отсчетов (оценка несостоятельная).

Автор: Бахурин Сергей www.dsplib.org

См. определение в файле psd_periodogram_cmplx.c строка 151

◆ psd_welch()

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
	)

Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом Уэлча.

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)

Функция рассчитывает спектральную плотность мощности $ X(f) $ выборки сигнала длительности $n $ отсчетов методом Уэлча:

\[ 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, \]

где $ w(m) $ – отсчёты оконной функции, $ F_s $ – частота дискретизации (Гц), $P = n/n_{\text{overlap}}$ – количество сегментов смещений выборки сигналов размера $n_{FFT}$,

$ U $ нормировочный коэффициент равный

\[ U = \sum_{m = 0}^{n-1} w^2(m), \]

Процедура разбиения исходной последовательности длительности n отсчетов на сегменты длины $n_{FFT}$ отсчетов, перекрывающихся с интервалом $n_{\text{overlap}}$ отсчетов, показан на следующем рисунке

Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого преобразования Фурье, для дискретной сетки частот от 0 Гц до $ F_s $ Гц (по умолчанию), или от $-F_s /2 $ до $F_s /2 $, если установлен флаг расчета двусторонней периодограммы.

Заметки: Периодограмма Уэлча возвращает смещенную, но состоятельную оценку СПМ.

Аргументы

[in]	x	Указатель на входной вектор комплексного сигнала $x(m)$, $ m = 0 \ldots n-1 $. Размер вектора `[n x 1]`.
[in]	n	Размер вектора входного сигнала. Также размер выходного вектора СПМ и вектора частоты также равны `n`.
[in]	win_type	Тип оконной функции, применяемой для модифицированной периодограммы. Подробнее смотри описание функции window.
[in]	win_param	Параметр оконной функциии. Данный параметр применяется только для парамтрических типов окон (смотри описание функции window). Для непараметрических функций игнорируется.
[in]	nfft	Размер перекрывающегося сегмента. Размер выходного вектора СПМ, и соответсвующего ей вектора частоты.
[in]	noverlap	Размер сдвига сегментов относительно друг друга (отсчетов). `noverlap = nfft` задает оценку без перекрытия сегментов. Обычно используют сдвиг равный половине размера сегмента `noverlap = nfft/2`.
[in]	pfft	Указатель на структуру fft_t. Указатель может быть `NULL`. В этом случае объект структуры будет создан внутри функции и удален перед завершением. Если предполагается многократный вызов функции, то рекомендуется создать объект fft_t и передавать в функцию, чтобы не создавать его каждый раз.
[in]	fs	частота дискретизации выборки исходного сигнала (Гц).
[in]	flag	Комбинация битовых флагов, задающих режим расчета: DSPL_FLAG_LOGMAG - СПМ считать в логарифмическом масштабе в единицах дБ/Гц DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2)
[in,out]	ppsd	Указатель на вектор СПМ рассчитанных по входному сигналу $x$. Размер вектора `[nfft x 1]`. Память должна быть выделена.
[in,out]	pfrq	Указатель на вектор частоты, соответствующей значениям рассчитанного вектора СПМ. Размер вектора `[nfft x 1]`. Указатель может быть `NULL`,в этом случае вектор частоты не рассчитывается и не возвращается.

Возвращает: RES_OK если расчет произведен успешно.
В противном случае код ошибки.

Пример периодограммных оценок СПМ для различной длины выборки сигнала:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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_welch(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_8192.txt");
    psd_plot(argc, argv,     "img/psd_welch_8192.png", 
                             "dat/psd_welch_8192.txt");
 
 
    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");
 
    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");
 
    dspl_free(hdspl);      /* free dspl handle */
    if(x)
      free(x);
    if(psd)
      free(psd);
    if(frq)
      free(frq);
  return 0;
}
 
 

Программа производит расчет СПМ сигнала, состоящего из двух комплексных гармоник на фоне белого гауссова шума. Расчет ведется по выборке сигнала длины 8192 отсчета.

Рассчитанные СПМ выводятся на графики:

nfft = 8192, noverlap = 4096:

nfft = 1024, noverlap = 512:

nfft = 256, noverlap = 128:

Можно видеть, что уменьшение nfft при фиксированной длительности сигнала позволяет уменьшить флуктуации шума и делает оценку состоятельной.

Автор: Бахурин Сергей www.dsplib.org

См. определение в файле psd_welch.c строка 160

◆ psd_welch_cmplx()

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
	)

Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом Уэлча.

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)

Функция рассчитывает спектральную плотность мощности $ X(f) $ выборки сигнала длительности $n $ отсчетов методом Уэлча:

\[ 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, \]

где $ w(m) $ – отсчёты оконной функции, $ F_s $ – частота дискретизации (Гц), $P = n/n_{\text{overlap}}$ – количество сегментов смещений выборки сигналов размера $n_{FFT}$,

$ U $ нормировочный коэффициент равный

\[ U = \sum_{m = 0}^{n-1} w^2(m), \]

Процедура разбиения исходной последовательности длительности n отсчетов на сегменты длины $n_{FFT}$ отсчетов, перекрывающихся с интервалом $n_{\text{overlap}}$ отсчетов, показан на следующем рисунке

Расчет спектральной плотности мощности ведется при помощи алгоритмов быстрого преобразования Фурье, для дискретной сетки частот от 0 Гц до $ F_s $ Гц (по умолчанию), или от $-F_s /2 $ до $F_s /2 $, если установлен флаг расчета двусторонней периодограммы.

Заметки: Периодограмма Уэлча возвращает смещенную, но состоятельную оценку СПМ.

Аргументы

[in]	x	Указатель на входной вектор комплексного сигнала $x(m)$, $ m = 0 \ldots n-1 $. Размер вектора `[n x 1]`.
[in]	n	Размер вектора входного сигнала. Также размер выходного вектора СПМ и вектора частоты также равны `n`.
[in]	win_type	Тип оконной функции, применяемой для модифицированной периодограммы. Подробнее смотри описание функции window.
[in]	win_param	Параметр оконной функциии. Данный параметр применяется только для парамтрических типов окон (смотри описание функции window). Для непараметрических функций игнорируется.
[in]	nfft	Размер перекрывающегося сегмента. Размер выходного вектора СПМ, и соответсвующего ей вектора частоты.
[in]	noverlap	Размер сдвига сегментов относительно друг друга (отсчетов). `noverlap = nfft` задает оценку без перекрытия сегментов. Обычно используют сдвиг равный половине размера сегмента `noverlap = nfft/2`.
[in]	pfft	Указатель на структуру fft_t. Указатель может быть `NULL`. В этом случае объект структуры будет создан внутри функции и удален перед завершением. Если предполагается многократный вызов функции, то рекомендуется создать объект fft_t и передавать в функцию, чтобы не создавать его каждый раз.
[in]	fs	частота дискретизации выборки исходного сигнала (Гц).
[in]	flag	Комбинация битовых флагов, задающих режим расчета: DSPL_FLAG_LOGMAG - СПМ считать в логарифмическом масштабе в единицах дБ/Гц DSPL_FLAG_PSD_TWOSIDED - двусторонняя СПМ (от -Fs/2 до Fs/2)
[in,out]	ppsd	Указатель на вектор СПМ рассчитанных по входному сигналу $x$. Размер вектора `[nfft x 1]`. Память должна быть выделена.
[in,out]	pfrq	Указатель на вектор частоты, соответствующей значениям рассчитанного вектора СПМ. Размер вектора `[nfft x 1]`. Указатель может быть `NULL`,в этом случае вектор частоты не рассчитывается и не возвращается.

Возвращает: RES_OK если расчет произведен успешно.
В противном случае код ошибки.

Пример периодограммных оценок СПМ для различной длины выборки сигнала:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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;
}
 
 

Программа производит расчет СПМ сигнала, состоящего из двух комплексных гармоник на фоне белого гауссова шума. Расчет ведется по выборке сигнала длины 8192 отсчета.

Рассчитанные СПМ выводятся на графики:

nfft = 8192, noverlap = 4096:

nfft = 1024, noverlap = 512:

nfft = 256, noverlap = 128:

Можно видеть, что уменьшение nfft при фиксированной длительности сигнала позволяет уменьшить флуктуации шума и делает оценку состоятельной.

Автор: Бахурин Сергей www.dsplib.org

См. определение в файле psd_welch_cmplx.c строка 159