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

Сгенерируйте двухканальный синусоидальный сигнал, дискретизированный с частотой дискретизации 1 кГц около 0,5 с. Синусоидальные частоты составляют 175 Гц и 400 Гц. Встройте сигнал в белый шум с отклонением $${\sigma }^2=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 года. Когда погода позволяла, Галилей регистрировал местоположение спутников. Используйте его наблюдения, чтобы оценить орбитальный период одного из спутников, Callisto.

Угловое положение Каллисто измеряется в минутах дуги. Отсутствующие данные из-за облачных условий задаются с помощью 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.5{day}^{-1}$$. Задайте коэффициент избыточной дискретизации 10. Выберите стандартную нормализацию Lomb-Scargle.

oneday = seconds(days(1));

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

f = f*oneday;

Периодограмма имеет один ясный максимум. Назовите пиковую частоту $$f_0$$. Постройте график периодограммы и аннотируйте пик.

[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=1024$$ выборки белого шума с отклонением $$\sigma =1$$, при скорости дискретизации 1 Гц. Вычислите степень спектр белого шума. Выберите нормализацию Lomb-Scargle и задайте коэффициент избыточной дискретизации $$ofac=15$$. Сбросьте генератор случайных чисел для повторяемых результатов.

rng default

N = 1024;
t = (1:N)';
wgn = randn(N,1);

ofac = 15;
[pwgn,f] = plomb(wgn,t,[],ofac,'normalized');

Проверьте, что оценка степени белого шума Ломба-Скаргля имеет экспоненциальное распределение с единичным средним. Постройте гистограмму значений pwgn и гистограмму множества экспоненциально распределенных случайных чисел, сгенерированных с помощью функции распределения $$f(z|1)=e^{-z}$$. Чтобы нормализовать гистограммы, напомните, что общее количество выборок периодограммы составляет $$N\times ofac/2$$. Задайте ширину интервала 0,25. Наложите график функции теоретического распределения.

dx = 0.25;
br = 0:dx:5;

Nf = N*ofac/2;

hpwgn = histcounts(pwgn,br)';

hRand = histcounts(-log(rand(Nf,1)),br)';

bend = br(1:end-1);

bar(bend,[hpwgn hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('wgn','Empirical pdf','Theoretical pdf')
hold off

Встройте в шум синусоидальный сигнал частоты 0,1 Гц. Используйте отношение сигнал/шум $$\xi =0.01$$. Задайте синусоидальную амплитуду, $$x_0$$, с использованием отношения $$x_0=\sigma \sqrt{2\xi }$$. Вычислите степень спектр сигнала и постройте его гистограмму вместе с эмпирическими и теоретическое распределение функциями.

SNR = 0.01;
x0 = sqrt(2*SNR);
sigsmall = wgn + x0*sin(2*pi/10*t);

[psigsmall,f] = plomb(sigsmall,t,[],ofac,'normalized');

hpsigsmall = histcounts(psigsmall,br)';

bar(bend,[hpsigsmall hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('sigsmall','Empirical pdf','Theoretical pdf')
hold off

Повторите расчет используя $$\xi =1$$. Теперь распределение заметно отличается.

SNR = 1;
x0 = sqrt(2*SNR);
siglarge = wgn + x0*sin(2*pi/10*t);

[psiglarge,f] = plomb(siglarge,t,[],ofac,'normalized');

hpsiglarge = histcounts(psiglarge,br)';

bar(bend,[hpsiglarge hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('siglarge','Empirical pdf','Theoretical pdf')
hold off

Сгенерируйте 100 выборки синусоидального сигнала со скоростью дискретизации 1 Гц. Задайте амплитуду 0,75 и частоту $$0.6/2\pi \approx 0.096Hz$$. Встройте сигнал в белый шум отклонения 0.902. Сбросьте генератор случайных чисел для повторяемых результатов.

rng default

X0 = 0.75;
f0 = 0.6;
vr = 0.902;

Nsamp = 100;
t = 1:Nsamp;
X = X0*cos(f0*(1:Nsamp))+randn(1,Nsamp)*sqrt(vr);

Отбрасывайте 10% выборок случайным образом. Постройте график сигнала.

X(randperm(Nsamp,Nsamp/10)) = NaN;

plot(t,X,'o')

Вычислите и постройте график нормированного спектра степени. Аннотировать уровни, которые соответствуют вероятностям ложного предупреждения 50%, 10%, 1% и 0,01%. Если вы генерируете много 90-выборочных сигналов белого шума с отклонением 0,902, то половина из них имеет один или несколько peaks выше 50% линии, 10% имеют один или несколько peaks выше 10% линии и так далее.

Pfa = [50 10 1 0.01]/100;
Pd = 1-Pfa;

[pxx,f,pth] = plomb(X,1:Nsamp,'normalized','Pd',Pd);

plot(f,pxx,f,pth*ones(size(f')))
xlabel('f')
text(0.3*[1 1 1 1],pth-.5,[repmat('P_{fa} = ',[4 1]) num2str(Pfa')])

В этом случае пик достаточно высок, что только около 0,01% возможных сигналов могут достичь его.

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

plomb(X,1:Nsamp,'normalized','Pd',Pd)

Когда задан вектор данных как единственный вход, plomb оценивает спектральную плотность степени с использованием нормированных частот.

Сгенерируйте 128 выборки синусоиды нормированной частоты $$\pi /2$$ рад/образец, встроенный в белый Гауссов шум отклонения 1/100.

t = (0:127)';
x = sin(2*pi*t/4);
x = x + randn(size(x))/10;

Оцените PSD с помощью процедуры Lomb-Scargle. Повторите расчет с periodogram.

[p,f] = plomb(x);
[pper,fper] = periodogram(x);

Постройте графики оценок PSD в децибелах. Проверьте, что результаты эквивалентны.

plot(f/pi,pow2db(p))
hold on
plot(fper/pi,pow2db(pper))

axis([0 1 -40 20])
xlabel('\omega / \pi')
ylabel('PSD')
legend('plomb','periodogram')

Оцените Lomb-Scargle PSD трехканального сигнала, состоящего из синусоидов. Задайте частоты как $$2\pi /3$$ рад/образец, $$\pi /2$$ рад/образец, и $$2\pi /5$$ рад/образец. Добавьте белый Гауссов шум отклонения 1/100. Использование plomb без выходных аргументов для вычисления и построения графика оценки PSD в децибелах.

x3 = [sin(2*pi*t/3) sin(2*pi*t/4) sin(2*pi*t/5)];
x3 = x3 + randn(size(x3))/10;

figure
plomb(x3)

Снова вычислите оценку PSD, но теперь удалите 25% данных случайным образом.

x3(randperm(numel(x3),0.25*numel(x3))) = NaN;

plomb(x3)

Если у вас нет временного вектора, используйте NaN's, чтобы задать отсутствующие выборки в сигнале.

Сгенерируйте 1024 выборки синусоиды нормированной частоты $$2\pi /5$$ рад/образец, встроенный в белый шум отклонения 1/100. Оцените спектральную плотность степени с помощью процедуры Ломба-Скаргля. Использование plomb без выходных аргументов для построения графика оценки.

t = (0:1023)';
x = sin(2*pi*t/5);
x = x + randn(size(x))/10;

[pxx,f] = plomb(x);

plomb(x)

Удалите каждую другую выборку путем назначения NaN's. Использование plomb для вычисления и построения графика PSD. Периодограмма пикует с той же частотой, потому что ось времени неизменна.

xnew = x;
xnew(2:2:end) = NaN;

plomb(xnew)

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

xdec = x(1:2:end);

plomb(xdec)