libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
Алгоритмы цифрового спектрального анализа

Функции

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_{\text{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_{\text{FFT}}\).

При использовании \(n_{\text{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.c строка 149

◆ 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 если расчет произведен успешно.
В противном случае код ошибки.

Пример оценок СПМ методом Бартлетта:
#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.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;
}

Программа производит расчет СПМ сигнала, состоящего из двух комплексных экспонент на фоне белого гауссова шума. Расчет ведется по выборк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.c строка 367

◆ 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.c строка 590

◆ 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.c строка 818

◆ 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.c строка 1060

◆ 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.c строка 1321

Структура параметров датчиков псевдослучайных чисел.
Definition: dspl.h:289
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:359
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)
Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом Уэлча.
Definition: psd.c:1321
void gnuplot_close(void *h)
Закрыть хэндл GNUPLOT.
Definition: gnuplot.c:325
int randn(double *x, int n, double mu, double sigma, random_t *prnd)
Генерация вектора нормально распределенных псевдослучайных чисел.
Definition: randgen.c:589
int psd_bartlett_cmplx(complex_t *x, int n, int nfft, fft_t *pfft, double fs, int flag, double *ppsd, double *pfrq)
Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом Бартлетта...
Definition: psd.c:367
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)
Непараметрическая оценка спектральной плотности мощности (СПМ) комплексного сигнала методом модифицир...
Definition: psd.c:818
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)
Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом Уэлча.
Definition: psd.c:1060
int psd_bartlett(double *x, int n, int nfft, fft_t *pfft, double fs, int flag, double *ppsd, double *pfrq)
Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом Бартлетт...
Definition: psd.c:149
void * dspl_load()
Произвести динамическую линковку и загрузить функции libdspl-2.0.
void gnuplot_cmd(void *h, char *cmd)
Функция посылает команду cmd пакету GNUPLOT, для построения или оформления графика,...
Definition: gnuplot.c:390
int writetxt(double *x, double *y, int n, char *fn)
Сохранить вещественные данные в текстовый файл
Definition: inout.c:491
int random_init(random_t *prnd, int type, void *seed)
Инициализация датчиков псевдослучайных чисел.
Definition: randgen.c:121
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:86
void dspl_free(void *handle)
Очищает связанную ранее динамическую библиотеку DSPL-2.0.
int psd_periodogram(double *x, int n, int win_type, double win_param, fft_t *pfft, double fs, int flag, double *ppsd, double *pfrq)
Непараметрическая оценка спектральной плотности мощности (СПМ) вещественного сигнала методом модифици...
Definition: psd.c:590
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:417
int gnuplot_create(int argc, char *argv[], int w, int h, char *fn_png, void **hplot)
Создать график GNUPLOT.
Definition: gnuplot.c:200