Метод 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 года. Когда погода позволила, Галилео записал местоположения спутников. Используйте его наблюдения, чтобы оценить орбитальный период одного из спутников, Каллисто.

Угловое положение Каллисто измеряется в минутах дуги. Недостающие данные из-за облачных условий заданы с помощью 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');

Проверьте, что оценка спектра мощности Lomb-Scargle белого шума имеет экспоненциальное распределение с модульным средним значением. Постройте гистограмму значений 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задавать недостающие выборки в сигнале.

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

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

[pxx,f] = plomb(x);

plomb(x)

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

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

plomb(xnew)

Удалите любую выборку путем субдискретизации. Функция теперь оценивает периодичность на дважды исходной частоте. Это - вероятно, не результат, который вы хотите.

xdec = x(1:2:end);

plomb(xdec)