Этот пример показывает, как оценить значимость синусоидального компонента в белом шуме с помощью 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-статистики Фишера. Следующий код 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] Вихерт, София, Константинос Фокианос и Корбиниан Штриммер. «Идентификация периодически выраженных транскриптов в данных временных рядов микрочипов». Биоинформатика. Том 20, 2004, стр. 5-20.