Тестирование значения на периодический компонент

В этом примере показано, как оценить значение синусоидального компонента в белом шуме с помощью g-статистической-величины Фишера. G-статистическая-величина Фишера является отношением самого большого значения периодограммы к сумме всех значений периодограммы по 1/2 интервала частоты, (0, Fs/2). Подробное описание g-статистической-величины и точного распределения может быть найдено в ссылках.

Создайте сигнал, состоящий из синусоиды на 100 Гц в белом Гауссовом шуме с нулевым средним значением и отклонением 1. Амплитуда синусоиды 0.25. Частота дискретизации составляет 1 кГц. Установите генератор случайных чисел на настройки по умолчанию для восстанавливаемых результатов.

rng default

Fs = 1e3;
t = 0:1/Fs:1-1/Fs;
x = 0.25*cos(2*pi*100*t)+randn(size(t));

Получите периодограмму сигнала с помощью periodogram. Исключите 0 и частота Найквиста (Fs/2). Постройте периодограмму.

[Pxx,F] = periodogram(x,rectwin(length(x)),length(x),Fs);
Pxx = Pxx(2:length(x)/2);

periodogram(x,rectwin(length(x)),length(x),Fs)

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate contains an object of type line.

Найдите максимальное значение периодограммы. G-статистическая-величина Фишера является отношением максимального значения периодограммы к сумме всех значений периодограммы.

[maxval,index] = max(Pxx);
fisher_g = Pxx(index)/sum(Pxx)
fisher_g = 0.0381

Максимальное значение периодограммы происходит на уровне 100 Гц, которые можно проверить путем нахождения частоты, соответствующей индексу максимального значения периодограммы.

F = F(2:end-1);
F(index)
ans = 100

Используйте дистрибутивные результаты, детализированные в ссылках, чтобы определить уровень значения, pval, из g-статистической-величины Фишера. Следующий код MATLAB® реализует уравнение (6) [2]. Используйте логарифм гамма функции, чтобы избежать переполнения когда вычислительные биномиальные коэффициенты.

N = length(Pxx);
nn = 1:floor(1/fisher_g);

I = (-1).^(nn-1).*exp(gammaln(N+1)-gammaln(nn+1)-gammaln(N-nn+1)).*(1-nn*fisher_g).^(N-1);

pval = sum(I)
pval = 2.0163e-06

P-значение меньше 0.00001, который указывает на значительный периодический компонент на уровне 100 Гц. Интерпретация g-статистической-величины Фишера осложнена присутствием других периодичностей. См. [1] для модификации, когда несколько периодичностей смогут присутствовать.

Ссылки

[1] Персиваль, Дональд Б. и Эндрю Т. Уолден. Спектральный анализ для физических приложений. Кембридж, Великобритания: Издательство Кембриджского университета, 1993.

[2] Wichert, София, Константинос Фокиэнос и Корбиниэн Стриммер. "Идентифицируя Периодически Описываемые Расшифровки стенограммы в Данных временных рядов Микромассивов". Биоинформатика. Издание 20, 2004, стр 5-20.

Смотрите также

|