このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

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 は単調増加しなければなりませんが、その間隔は均一である必要はありません。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 を利用したゼロ パディング手法に類似しています。ここでも、pxx が round(fmax/f_min) 周波数点で返されますが、この場合に考慮される最小周波数は 1/(ofac × N × t_s) となります。ofac の既定値は 4 です。

[pxx,fvec] = plomb(___,fvec) は fvec で指定されている周波数で x の PSD を推定します。fvec は少なくとも 2 つの要素をもたなければなりません。2 番目の出力引数は入力 fvec と同じです。

この構文を使用している場合、最大周波数またはオーバーサンプリング係数を指定することはできません。

[___] = plomb(___,spectrumtype) では、ピリオドグラムの正規化を指定します。

• spectrumtype を 'psd' に設定するか、何も指定せずに、パワー スペクトル密度として pxx を求めます。

• spectrumtype を 'power' に設定して入力信号のパワー スペクトルを求めます。

• spectrumtype を 'normalized' に設定して、x の分散を 2 倍してスケーリングされる、標準の Lomb-Scargle ピリオドグラムを得ます。

[___,pth] = plomb(___,'Pd',pdvec) はパワーレベルのしきい値 pth を返します。pth よりも大きな値をもつピークは pdvec の確率で、ランダムな変動の結果ではなく真の信号ピークです。pdvec はベクトルにすることができます。pdvec の各要素は 0 より大きく 1 未満でなければなりません。pth の各行は、pdvec の各要素に対応しています。pth は x と同じチャネル数をもちます。このオプションは出力周波数を fvec で指定した場合には使用できません。

[pxx,w] = plomb(x) では等間隔に並んだ一連の正規化周波数 w で評価された x の PSD 推定を返し、ナイキスト区間をカバーします。NaN を使用して欠損サンプルを指定します。上記のオプションはすべて、正規化周波数で使用できます。オプションを使用するには、2 番目の入力として空配列を指定します。

出力引数を設定せずに plomb(___) を使用すると、現在の Figure ウィンドウに Lomb-Scargle ピリオドグラム PSD 推定がプロットされます。

例

不定間隔でサンプリングされた信号と欠損サンプルのある信号

Lomb-Scargle 法では、不定間隔でサンプリングされた信号と欠損サンプルのある信号を処理できます。

1 kHz で約 0.5 秒間サンプリングされた 2 チャネルの正弦波信号を生成します。正弦波の周波数は 175 Hz および 400 Hz です。分散が $${\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')

元の信号をサンプリングしますが、ジッター (不確かさ) を時間測定値に追加することで、サンプリングを不等間隔に行います。最初の時点は、ゼロになるまで継続します。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 つの衛星の動きを観察しました。ガリレオは天候が許す限り衛星の位置を記録しました。この観察記録を使用して衛星の 1 つであるカリストの軌道周期を推定します。

カリストの角度位置は分角で測定されています。曇天などの気象条件によって記録できなかったデータは、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 つの衛星の中で最も外側にあるカリストの軌道周期を求めます。

ガリレオはカリストの角度位置を分角で測定しています。曇天などの気象条件による空白箇所がいくつかあります。観測時間の 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;

ピリオドグラムには明確な最大値が 1 つあります。そのピーク周波数を $$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)

不規則なサンプリングとエイリアシング

1 Hz で 100 秒間 0.8 Hz の正弦波をサンプリングします。分散 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 Hz までの周波数を含めます。オーバーサンプリング係数を 2 に指定します。PSD をプロットします。

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

ofac = 2;

plomb(n,tn,1,ofac)

エイリアシングが消えています。不規則なサンプリングでは、時間間隔を縮小することで有効なサンプル レートが増加します。

0.8 Hz 前後の周波数を拡大します。間隔が 0.001 Hz の細かいグリッドを使用します。この場合、オーバーサンプリング係数または最大周波数を指定することはできません。

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

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

指数分布

1 Hz のサンプル レートにおいて、分散が $$\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$$ であることを思い出してください。ビンの幅を 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 Hz の正弦波信号をノイズに組み込みます。$$\xi =0.01$$ の S/N 比を使用します。正弦波の振幅 $$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

誤警報率

サンプル レートが 1 Hz の正弦波信号のサンプルを 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% のラインより高いピークを 1 つ以上もつ、10% は、10% のラインより高いピークを 1 つ以上もつ、などです。

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$$ ラジアン/サンプルの正弦波の 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 チャネルの信号の PSD を Lomb-Scargle 法により推定します。周波数を $$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/100 のホワイト ノイズに組み込まれる、正規化周波数 $$2\pi /5$$ ラジアン/サンプルの正弦波の 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 を代入して、1 つおきにサンプルを取り除きます。plomb を使用して PSD を計算してプロットします。ピリオドグラムは、時間軸が変わらないため、同じ周波数でピークになります。

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

plomb(xnew)

ダウンサンプリングで 1 つおきにサンプルを取り除きます。ここで関数は元の周波数の 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 は単調増加しなければなりませんが、その間隔は均一である必要はありません。t には NaN または NaT を含めることができます。これらの値は欠損データとして取り扱われ、スペクトルの計算からは除外されます。

データ型: single | double | datetime | duration

fs — サンプル レート

正のスカラー

サンプル レート。正のスカラーで指定します。サンプル レートは単位時間あたりのサンプル数です。時間の単位が秒の場合、サンプル レートの単位は Hz です。

データ型: 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' を指定すると、x の分散を 2 倍して pxx をスケーリングします。

pdvec — 検出の確率

スカラー | ベクトル

検出の確率。'Pd' とスカラーまたは 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]オフセットは、測定時点にのみ依存し、その時点が等間隔に配置されたときにゼロになります。

入力信号がホワイト ガウス ノイズから構成されている場合、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++ コードを生成します。

使用上の注意および制限:

コード生成では、timetable はサポートされていません。

バージョン履歴

R2014b で導入

参考

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

トピック

• 欠損サンプルのある信号の周期性の検出
• 不等間隔サンプル信号のスペクトル解析