Метод Ломба-Скаргла может обрабатывать сигналы, которые были отобраны неравномерно или имеют отсутствующие выборки.

Генерируют двухканальный синусоидальный сигнал, дискретизированный на частоте 1 кГц в течение примерно 0,5 С. Синусоидальные частоты составляют 175 Гц и 400 Гц. Встраивайте сигнал в белый шум с дисперсией $${}^{\mathrm{start2}}\mathrm{=}\mathrm{}$$1/4.

Fs = 1000;
f0 = 175;
f1 = 400;

t = 0:1/Fs:0.5;

wgn = randn(length(t),2)/2;

sigOrig = sin(2*pi*[f0;f1]*t)' + wgn;

Вычислите и постройте график периодограммы сигнала. Использовать periodogram с настройками по умолчанию.

periodogram(sigOrig,[],[],Fs)

axisLim = axis;
title('Periodogram')

Использовать plomb с настройками по умолчанию для оценки и построения графика PSD сигнала. Используйте границы оси из предыдущего графика.

plomb(sigOrig,t)

axis(axisLim)
title('Lomb-Scargle')

Предположим, что в сигнале отсутствует 10% исходных выборок. Место NaNs в случайных местах для моделирования отсутствующих точек данных. Использовать plomb для оценки и построения графика PSD сигнала с отсутствующими выборками.

sinMiss = sigOrig;

misfrac = 0.1;
nTime = length(t)*2;

sinMiss(randperm(nTime,round(nTime*misfrac))) = NaN;

plomb(sinMiss,t)

axis(axisLim)
title('Missing Samples')

Выполните выборку исходного сигнала, но сделайте выборку неравномерной, добавив дрожание (неопределенность) к измерениям времени. Первое мгновение продолжает оставаться на нуле. Использовать plomb для оценки и построения графика PSD неравномерно дискретизированного сигнала.

tirr = t + (1/2-rand(size(t)))/Fs/2;
tirr(1) = 0;

sinIrreg = sin(2*pi*[f0;f1]*tirr)' + wgn;

plomb(sinIrreg,tirr)

axis(axisLim)
title('Nonuniform Sampling')

Галилей Галилей наблюдал движение четырех крупнейших спутников Юпитера зимой 1610 года. Когда погода позволила, Галилей зафиксировал местонахождение спутников. Используйте его наблюдения, чтобы оценить орбитальный период одного из спутников, Каллисто.

Угловое положение Каллисто измеряют в минутах дуги. Недостающие данные из-за облачности уточняются с помощью NaNs. Первое наблюдение датировано 15 января. Создать datetime массив времен наблюдения.

yg = [10.5 NaN 11.5 10.5 NaN NaN NaN -5.5 -10.0 -12.0 -11.5 -12.0 -7.5 ...
    NaN NaN NaN NaN 8.5 12.5 12.5 10.5 NaN NaN NaN -6.0 -11.5 -12.5 ...
    -12.5 -10.5 -6.5 NaN 2.0 8.5 10.5 NaN 13.5 NaN 10.5 NaN NaN NaN ...
    -8.5 -10.5 -10.5 -10.0 -8.0]';

obsv = datetime(1610,1,14+(1:length(yg)));

plot(yg,'o')

ax = gca;
nights = [1 18 32 46];
ax.XTick = nights;
ax.XTickLabel = char(obsv(nights));
grid

Оценка спектра мощности данных с помощью plomb. Укажите коэффициент избыточной выборки, равный 10. Выражайте результирующие частоты в обратные дни.

[pxx,f] = plomb(yg,obsv,[],10,'power');
f = f*86400;

Использовать findpeaks для определения местоположения единственного выдающегося пика спектра. Постройте график спектра мощности и отобразите пик.

[pk,f0] = findpeaks(pxx,f,'MinPeakHeight',10);

plot(f,pxx,f0,pk,'o')
xlabel('Frequency (day^{-1})')
title('Power Spectrum and Prominent Peak')
grid

Определите орбитальный период Каллисто (в сутках) как обратный частоте максимальной энергии. Результат отличается менее чем на 1% от значения, опубликованного НАСА.

Period = 1/f0

Period = 16.6454

NASA = 16.6890184;
PercentDiscrep = (Period-NASA)/NASA*100

PercentDiscrep = -0.2613

Галилео Галилей обнаружил четыре крупнейших спутника Юпитера в январе 1610 года и записывал их местоположение каждую ясную ночь до марта того же года. Используйте данные Галилея, чтобы найти орбитальный период Каллисто, самого дальнего из четырех спутников.

Наблюдения Галилея за угловым положением Каллисто в минутах дуги. Есть несколько разрывов из-за мутных условий. Создать duration массив времен наблюдения.

t = [0 2 3 7 8 9 10 11 12 17 18 19 20 24 25 26 27 28 29 31 32 33 35 37 ...
    41 42 43 44 45]';
td = days(t);

yg = [10.5 11.5 10.5 -5.5 -10.0 -12.0 -11.5 -12.0 -7.5 8.5 12.5 12.5 ...
    10.5 -6.0 -11.5 -12.5 -12.5 -10.5 -6.5 2.0 8.5 10.5 13.5 10.5 -8.5 ...
    -10.5 -10.5 -10.0 -8.0]';

plot(td,yg,'o')

Использовать plomb для вычисления периодограммы данных. Оцените спектр мощности до частоты $$0\mathrm{.}{}^{\mathrm{}}$$5дня-1. Укажите коэффициент избыточной выборки, равный 10. Выберите стандартную нормализацию Lomb-Scargle.

oneday = seconds(days(1));

[pxx,f] = plomb(yg,td,0.5/oneday,10,'normalized');

f = f*oneday;

Периодограмма имеет один четкий максимум. Назовите пиковую частоту $${}_{\mathrm{f0}}$$. Постройте график периодограммы и аннотируйте пик.

[pmax,lmax] = max(pxx);
f0 = f(lmax);

plot(f,pxx,f0,pmax,'o')
xlabel('Frequency (day^{-1})')

Использование линейных наименьших квадратов для соответствия данным функции формы

Параметрами фитинга являются амплитуды A, B и C.

ft = 2*pi*f0*t;

ABC = [ones(size(ft)) cos(ft) sin(ft)] \ yg

ABC = 3×1

    0.4243
   10.4444
    6.6137

Используйте параметры фитинга для построения функции фитинга на 200-точечном интервале. Постройте график данных и наложите аппроксимацию.

x = linspace(t(1),t(end),200)';
fx = 2*pi*f0*x;

y = [ones(size(fx)) cos(fx) sin(fx)] * ABC;

plot(td,yg,'o',days(x),y)

Образец синусоиды 0,8 Гц при 1 Гц в течение 100 с. Встраивание синусоиды в белый шум с дисперсией 1/100. Сбросьте генератор случайных чисел для воспроизводимых результатов.

f0 = 0.8;

rng default
wgn = randn(1,100)/10;

ts = 1:100;
s = sin(2*pi*f0*ts) + wgn;

Вычислите и постройте график оценки спектральной плотности мощности до частоты дискретизации. Укажите коэффициент избыточной выборки, равный 10.

plomb(s,ts,1,10)

Существует наложение, потому что частота синусоиды больше частоты Найквиста.

Повторите расчет, но теперь случайным образом выполните выборку синусоиды. Включить частоты до 1 Гц. Укажите коэффициент избыточной выборки, равный 2. Постройте график PSD.

tn = sort(100*rand(1,100));
n = sin(2*pi*f0*tn) + wgn;

ofac = 2;

plomb(n,tn,1,ofac)

Псевдоним исчезает. Нерегулярная выборка увеличивает эффективную частоту выборки за счет сокращения некоторых временных интервалов.

Увеличьте изображение частот около 0,8 Гц. Используйте тонкую сетку с интервалом 0,001 Гц. В этом случае нельзя указать коэффициент избыточной дискретизации или максимальную частоту.

df = 0.001;
fvec = 0.7:df:0.9;

hold on
plomb(n,tn,fvec)
legend('ofac = 2','df = 0.001')

Генерация $$N\mathrm{=}\mathrm{}\mathrm{}\mathrm{}$$1024 выборок белого шума с $$\mathrm{дисперсией}\mathrm{\lambda }$$= 1, учитывая частоту выборок 1 Гц. Вычислите спектр мощности белого шума. Выберите нормализацию Lomb-Scargle и укажите