Проверка значимости для периодического компонента
Этот пример показывает, как оценить значимость синусоидального компонента в белом шуме с помощью 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)
Найдите максимальное значение периодограммы. 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-statistic Фишера. Следующий код 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] Персиваль, Дональд Б. и Эндрю Т. Уолден. Спектральный анализ для физических применений. Кембридж, Великобритания: Cambridge University Press, 1993.
[2] Wichert, Sofia, Konstantinos Fokianos, and Korbinian Strimmer. «Идентификация периодически выражаемых транскриптов в Микромассив Данных временных рядов». Биоинформатика. Том 20, 2004, стр. 5-20.