plomb

Lomb-Scargle 주기도

구문

[pxx,f] = plomb(x,t)

[pxx,f] = plomb(x,fs)

[pxx,f] = plomb(___,fmax)

[pxx,f] = plomb(___,fmax,ofac)

[pxx,fvec] = plomb(___,fvec)

[___] = plomb(___,spectrumtype)

[___,pth] = plomb(___,'Pd',pdvec)

[pxx,w] = plomb(x)

plomb(___)

설명

[pxx,f] = plomb(x,t)는 t에 지정된 시점에 샘플링된 신호 x의 Lomb-Scargle 파워 스펙트럼 밀도(PSD) 추정값 pxx를 반환합니다. t는 단조 증가(Monotonically Increasing)해야 하지만, 간격이 균일할 필요는 없습니다. t의 모든 요소는 음수가 아니어야 합니다. pxx는 f로 반환된 주파수에서 계산됩니다.

• x가 벡터이면 단일채널로 처리됩니다.

• x가 행렬이면 plomb는 각 열의 PSD를 개별적으로 계산하고 그 값을 pxx의 해당 열로 반환합니다.

x 또는 t에는 NaN 또는 NaT가 포함될 수 있습니다. 단, 이러한 값은 누락된 데이터로 처리되므로 스펙트럼 계산에서 제외됩니다.

[pxx,f] = plomb(x,fs)는 신호가 샘플 레이트 fs로 균일하게 샘플링되었지만 샘플이 누락된 경우를 처리합니다. NaN을 사용하여 누락된 데이터를 지정합니다.

[pxx,f] = plomb(___,fmax)는 위 구문의 입력 인수를 사용하여 최대 주파수 fmax까지 PSD를 추정합니다. 신호가 NaN이 아닌 N개의 시점에서 샘플링되고 Δt가 첫 번째 시점과 마지막 시점의 시간 차분인 경우, pxx는 round(fmax/f_min)개의 지점에서 반환됩니다. 여기서 f_min = 1/(4 × N × t_s)는 pxx가 계산된 가장 작은 주파수이고, 평균 샘플 시간은 t_s = Δt/(N – 1)입니다. fmax는 디폴트 값이 1/(2 × t_s)이며, 균일하게 샘플링된 신호인 경우 나이퀴스트 주파수에 해당합니다.

[pxx,f] = plomb(___,fmax,ofac)는 정수 오버샘플링 인자 ofac를 지정합니다. ofac를 사용하여 스펙트럼을 보간하거나 평활화하는 방식은 FFT 기반 방법의 0 채우기 기법과 유사합니다. 여기서도 pxx는 round(fmax/f_min)개의 주파수 지점에서 반환되지만, 이 경우에 고려되는 최소 주파수는 1/(ofac×N×t_s)입니다. ofac의 디폴트 값은 4입니다.

[pxx,fvec] = plomb(___,fvec)는 fvec에 지정된 주파수에서 x의 PSD를 추정합니다. fvec는 적어도 2개의 요소를 가져야 합니다. 두 번째 출력 인수는 입력 인수 fvec와 동일합니다.

이 구문을 사용할 경우, 최대 주파수나 오버샘플링 인자를 지정할 수 없습니다.

[___] = plomb(___,spectrumtype)은 주기도의 정규화를 지정합니다.

• pxx를 파워 스펙트럼 밀도로 구하려면 spectrumtype을 'psd'로 설정하거나 아무 값도 지정하지 마십시오.

• 입력 신호의 파워 스펙트럼을 구하려면 spectrumtype을 'power'로 설정하십시오.

• x의 분산의 2배로 스케일링된 표준 Lomb-Scargle 주기도를 구하려면, spectrumtype을 'normalized'로 설정하십시오.

[___,pth] = plomb(___,'Pd',pdvec)는 파워 수준 임계값 pth를 반환하며, 이 경우 pth보다 큰 값을 갖는 피크가 불규칙적인 변동으로 인한 결과가 아닌, 실제 신호 피크일 확률 pdvec를 갖습니다. pdvec는 벡터일 수 있습니다. pdvec의 모든 요소는 0보다 크고 1보다 작아야 합니다. pth의 각 행은 pdvec의 요소에 대응합니다. pth는 x와 채널 수가 동일합니다. 이 옵션은 fvec에 출력 주파수를 지정한 경우 사용할 수 없습니다.

[pxx,w] = plomb(x)는 나이퀴스트 구간에 걸쳐 있는 균일한 간격의 정규화 주파수 w에서 계산한 x의 PSD 추정값을 반환합니다. 누락된 샘플을 지정하려면 NaN을 사용하십시오. 위의 모든 옵션은 정규화 주파수에 사용 가능합니다. 이러한 옵션에 액세스하려면 두 번째 입력값으로 빈 행렬을 지정하십시오.

plomb(___)(출력 인수 없음)는 현재 Figure 창에 Lomb-Scargle 주기도 PSD 추정값을 플로팅합니다.

예제

불규칙적으로 샘플링된 신호 및 누락된 샘플이 있는 신호

Lomb-Scargle 방법은 균일하지 않게 샘플링되었거나 누락된 샘플이 있는 신호를 처리할 수 있습니다.

약 0.5초 동안 1kHz로 샘플링된 2채널 정현파 신호를 생성합니다. 정현파 주파수는 175Hz와 400Hz입니다. 이 신호에 분산이 $${\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%가 누락되었다고 가정하겠습니다. NaN을 임의의 위치에 배치하여 누락된 데이터 점을 시뮬레이션합니다. 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')

원래 신호를 샘플링하되 시간 측정값에 지터(불확실성)를 추가하여 불균일하게 샘플링합니다. 첫 번째 시점은 계속해서 0에 있습니다. 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년 겨울 동안 목성의 가장 큰 위성 4개의 움직임을 관측했습니다. 날씨가 허락하는 경우 갈릴레오는 위성의 위치를 기록했습니다. 갈릴레오의 관측을 토대로 위성 중 하나인 칼리스토의 궤도 주기를 추정합니다.

칼리스토의 각위치는 분각(minutes of arc) 단위로 측정되었습니다. 구름 낀 날씨로 인해 누락된 데이터는 NaN으로 지정되었습니다. 첫 번째 관측 날짜는 1월 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으로 지정합니다. 결과로 생성된 주파수를 1일 초수의 역수로 표현합니다.

[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

칼리스토의 궤도 주기(단위: 일)를 최대 에너지를 갖는 주파수의 역으로 구합니다. 그 결과값은 NASA에서 공개한 값과 1% 미만의 차이가 납니다.

Period = 1/f0

Period = 16.6454

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

PercentDiscrep = -0.2613

불규칙적으로 샘플링된 데이터 세트의 주기도

갈릴레이 갈릴레오는 1610년 1월에 목성의 가장 큰 위성 4개를 발견했고, 그 해 3월까지 밤하늘이 깨끗한 날이면 항상 그 위성들의 위치를 기록했습니다. 갈릴레오의 데이터를 사용하여 4개 위성 중 가장 바깥쪽에 있는 칼리스토의 궤도 주기를 구합니다.

갈릴레오는 칼리스토의 각위치를 분각(minutes of arc) 단위로 측정했습니다. 구름이 낀 날이 있었기 때문에 여러 개의 공백이 있습니다. 관측 시간값으로 구성된 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})')

선형 최소제곱을 사용하여 데이터에 다음 형식의 함수를 피팅합니다.

$$y(t)=A+B\cos 2\pi f_{0}t+C\sin 2\pi f_{0}t.$$

피팅 파라미터는 진폭 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)

불규칙적인 샘플링과 에일리어싱

100초 동안 1Hz로 0.8Hz 정현파를 샘플링합니다. 이 정현파에 분산이 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)

정현파 주파수가 나이퀴스트 주파수보다 크기 때문에 에일리어싱이 발생합니다.

계산을 반복하되 이번에는 임의 시간에 정현파를 샘플링합니다. 1Hz 이하의 주파수만 포함합니다. 오버샘플링 인자를 2로 지정합니다. PSD를 플로팅합니다.

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

ofac = 2;

plomb(n,tn,1,ofac)

에일리어싱이 사라집니다. 불규칙적인 샘플링으로 일부 시간 구간이 축소되어 유효 샘플 레이트가 증가되었습니다.

0.8Hz 주변 주파수를 확대합니다. 간격이 0.001Hz인 조밀한 그리드를 사용합니다. 이 경우에는 오버샘플링 인자나 최대 주파수를 지정할 수 없습니다.

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

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

지수 분포

샘플 레이트 1Hz에서 분산이 $$\sigma =1$$인 백색 잡음의 $$N=1024$$개 샘플을 생성합니다. 백색 잡음의 파워 스펙트럼을 계산합니다. 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$$라는 점에 유의하십시오. Bin 너비를 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.1Hz인 정현파 신호에 포함시킵니다. 신호 대 잡음비로 $$\xi =0.01$$을 사용합니다. 관계식 $$x_0=\sigma \sqrt{2\xi }$$를 사용하여 정현파 진폭 $$x_0$$을 지정합니다. 신호의 파워 스펙트럼을 계산하고, 그 히스토그램을 경험적 분포 함수와 이론적 분포 함수와 함께 플로팅합니다.

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

오경보 확률

1Hz의 샘플 레이트에서 100개의 정현파 신호 샘플을 생성합니다. 진폭은 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%에 해당하는 레벨에 주석을 표시합니다. 분산 0.902를 갖는 90개의 많은 백색 잡음 신호를 생성하는 경우, 그 잡음의 절반은 50% 선을 넘는 피크를 한 개 이상 갖고 10%는 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)

잡음이 있는 정현파에 대한 Lomb-Scargle 주기도

데이터 벡터가 유일한 입력값으로 주어진 경우, plomb는 정규화 주파수를 사용하여 파워 스펙트럼 밀도를 추정합니다.

분산이 1/100인 백색 가우스 잡음이 포함되어 있는 정규화 주파수 $$\pi /2$$ rad/sample로 구성된 128개의 정현파 샘플을 생성합니다.

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

Lomb-Scargle 절차를 사용하여 PSD를 추정합니다. 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')

정현파로 구성된 3채널 신호에 대한 Lomb-Scargle PSD를 추정합니다. 주파수를 $$2\pi /3$$ rad/sample, $$\pi /2$$ rad/sample, $$2\pi /5$$ rad/sample로 지정합니다. 분산이 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/100인 백색 가우스 잡음이 포함되어 있는 정규화 주파수 $$2\pi /5$$ rad/sample로 구성된 1024개의 정현파 샘플을 생성합니다. 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을 할당하여 샘플 2개마다 하나씩 샘플을 제거합니다. plomb를 사용하여 PSD를 계산하고 플로팅합니다. 시간 축이 변경되지 않았기 때문에 주기도는 동일한 주파수에서 피크에 도달합니다.

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

plomb(xnew)

다운샘플링하여 샘플 2개마다 하나씩 샘플을 제거합니다. 이제 함수는 원래 주파수의 두 배에 해당하는 주파수에서 주기성을 추정합니다. 원하는 결과가 반환되지 않을 수 있습니다.

xdec = x(1:2:end);

plomb(xdec)

입력 인수

x — 입력 신호

벡터 | 행렬

입력 신호로, 벡터나 행렬로 지정됩니다. x가 벡터이면 단일채널로 처리됩니다. x가 행렬이면 plomb는 각 열의 PSD 추정값을 개별적으로 계산하고 그 값을 pxx의 해당 열로 반환합니다. x에는 NaN이 포함될 수 있습니다. NaN은 누락된 데이터로 처리되므로 스펙트럼 계산에서 제외됩니다.

데이터형: single | double

t — 시점

음이 아닌 실수형 벡터 | datetime형 배열 | duration형 배열

시점으로, 음이 아닌 실수형 벡터, datetime형 배열 또는 duration형 배열로 지정됩니다. t는 단조 증가(Monotonically Increasing)해야 하지만 간격이 균일할 필요는 없습니다. t는 NaN 또는 NaT를 포함할 수 있습니다. 단, 이러한 값은 누락된 데이터로 처리되므로 스펙트럼 계산에서 제외됩니다.

데이터형: single | double | datetime | duration

fs — 샘플 레이트

양의 스칼라

샘플 레이트로, 양의 스칼라로 지정됩니다. 샘플 레이트는 단위 시간당 샘플 개수입니다. 시간 단위가 초이면 샘플 레이트의 단위는 헤르츠입니다.

데이터형: single | double

fmax — 최대 주파수

양의 스칼라

최대 주파수로, 양의 스칼라로 지정됩니다. fmax는 나이퀴스트 주파수보다 높을 수 있습니다.

데이터형: single | double

ofac — 오버샘플링 인자

4 (디폴트 값) | 양의 정수 스칼라

오버샘플링 인자로, 양의 정수 스칼라로 지정됩니다.

데이터형: single | double

fvec — 입력 주파수

벡터

입력 주파수로, 벡터로 지정됩니다. fvec는 적어도 2개의 요소를 가져야 합니다.

데이터형: single | double

spectrumtype — 파워 스펙트럼 스케일링

'psd' (디폴트 값) | 'power' | 'normalized'

파워 스펙트럼 스케일링으로, 'psd', 'power' 또는 'normalized' 중 하나로 지정됩니다. spectrumtype을 생략하거나 'psd'를 지정하면 파워 스펙트럼 밀도 추정값이 반환됩니다. 'power'를 지정하면 PSD의 각 추정값이 윈도우의 등가 잡음 대역폭으로 스케일링됩니다. 각 주파수에서 전력의 추정값을 구하려면 'power'를 지정하십시오. 'normalized'를 지정하면 pxx가 x의 분산의 두 배로 스케일링됩니다.

pdvec — 검출 확률

스칼라 | 벡터

검출 확률로, 'Pd'와 함께 0과 1 사이(0과 1 불포함)의 실수 값으로 구성된 스칼라 또는 벡터가 쉼표로 구분되어 지정됩니다. 검출 확률은 스펙트럼의 피크가 불규칙적인 변동으로 인한 것이 아닐 확률입니다.

데이터형: single | double

출력 인수

pxx — Lomb-Scargle 주기도

벡터 | 행렬

Lomb-Scargle 주기도로, 벡터나 행렬로 반환됩니다. 입력 신호 x가 벡터이면 pxx는 벡터입니다. x가 행렬이면, 이 함수는 x의 각 열을 개별 채널로 처리하고 각 채널의 주기도를 계산합니다.

f — 주파수

벡터

주파수로, 벡터로 반환됩니다.

데이터형: single | double

w — 정규화 주파수

벡터

정규화 주파수로, 벡터로 반환됩니다.

데이터형: single | double

pth — 전력 수준 임계값

벡터 | 행렬

전력 수준 임계값으로, 벡터나 행렬로 반환됩니다. 전력 수준 임계값은 스펙트럼의 피크가 불규칙적인 변동으로 인한 것이라는 가정에서 배제(확률 pdvec를 가짐)되기 위해 스펙트럼의 피크가 초과해야 하는 진폭입니다. pth의 각 행은 pdvec의 요소에 대응합니다. pth는 x와 채널 수가 동일합니다.

데이터형: single | double

세부 정보

Lomb-Scargle 주기도

Lomb-Scargle 주기도를 사용하면 불규칙하고 불균등하게 샘플링된 랜덤 데이터에서 약한 주기 신호를 찾고 테스트할 수 있습니다.

시간 t_k에서 가져온 N개의 관측값 x_k를 살펴보겠습니다. 여기서 k = 1, …, N입니다. Lomb-Scargle 주기도는 다음으로 정의됩니다.[2]

$$P_{\text{LS}}(f)=\frac{1}{2{\sigma }^2}\left\{\frac{{\left[{\displaystyle {\sum }_{k=1}^N(x_k-\overline{x})\cos \left(2\pi f(t_k-\tau )\right)}\right]}^2}{{\displaystyle {\sum }_{k=1}^N{\cos }^2\left(2\pi f(t_k-\tau )\right)}}+\frac{{\left[{\displaystyle {\sum }_{k=1}^N(x_k-\overline{x})\sin \left(2\pi f(t_k-\tau )\right)}\right]}^2}{{\displaystyle {\sum }_{k=1}^N{\sin }^2\left(2\pi f(t_k-\tau )\right)}}\right\},$$

여기서

$$\overline{x}=\frac{1}{N}{\displaystyle {\sum }_{k=1}^N{x}_k}$$

및

$${\sigma }^2=\frac{1}{N-1}{\displaystyle {\sum }_{k=1}^N{(x_k-\overline{x})}^2}$$

는 각각 데이터의 평균과 분산입니다.

시간 오프셋 τ는 다음이 성립되도록 선택됩니다.

$$\tan \left(2(2\pi f)\tau \right)=\frac{{\displaystyle {\sum }_{k=1}^N\sin \left(2(2\pi f)t_k\right)}}{{\displaystyle {\sum }_{k=1}^N\cos \left(2(2\pi f)t_k\right)}}$$

이는 계산된 스펙트럼의 시간 불변을 보장하기 위함입니다. 시간 측정값에 변위 t_k → t_k + T가 발생하면 오프셋에서도 다음과 같이 동일한 변위 τ → τ + T가 발생합니다. 또한 오프셋을 이와 같이 선택하면, "데이터에 대한 사인파 피팅 잔차의 제곱합이 최소화되는 주파수에서도 주기도의 최댓값이 발생 가능"합니다. [4] 오프셋은 측정 시간의 영향만 받으며 시간 간격이 균일할 경우 0이 됩니다.

입력 신호가 백색 가우스 잡음으로 구성된 경우 P_LS(f)는 단위 평균을 갖는 지수 확률 분포를 따릅니다[3].

참고 문헌

[1] Horne, James H., and Sallie L. Baliunas. "A Prescription for Period Analysis of Unevenly Sampled Time Series." Astrophysical Journal. Vol. 302, 1986, pp. 757–763.

[2] Lomb, Nicholas R. "Least-Squares Frequency Analysis of Unequally Spaced Data." Astrophysics and Space Science. Vol. 39, 1976, pp. 447–462.

[3] Press, William H., and George B. Rybicki. "Fast Algorithm for Spectral Analysis of Unevenly Sampled Data." Astrophysical Journal. Vol. 338, 1989, pp. 277–280.

[4] Scargle, Jeffrey D. "Studies in Astronomical Time Series Analysis. II. Statistical Aspects of Spectral Analysis of Unevenly Spaced Data." Astrophysical Journal. Vol. 263, 1982, pp. 835–853.

확장 기능

C/C++ 코드 생성

MATLAB® Coder™를 사용하여 C 코드나 C++ 코드를 생성할 수 있습니다.

사용법 관련 참고 및 제한 사항:

코드 생성에는 타임테이블이 지원되지 않습니다.

버전 내역

R2014b에 개발됨

참고 항목

bandpower | pburg | pcov | peig | periodogram | pmcov | pmtm | pmusic | pwelch | pyulear | spectrogram

도움말 항목

• Detect Periodicity in a Signal with Missing Samples
• 불균일하게 샘플링된 신호의 스펙트럼 분석

×

예제 열기

이 예제의 수정된 버전이 있습니다. 사용자가 편집한 내용을 반영하여 이 예제를 여시겠습니까?

아니요, 수정된 버전을 덮어씁니다. 예