pspectrum

周波数領域および時間-周波数領域内の信号の解析

構文

p = pspectrum(x)

p = pspectrum(x,fs)

p = pspectrum(x,t)

p = pspectrum(___,type)

p = pspectrum(___,Name,Value)

[p,f] = pspectrum(___)

[p,f,t] = pspectrum(___,'spectrogram')

[p,f,pwr] = pspectrum(___,'persistence')

pspectrum(___)

説明

p = pspectrum(x) は、x のパワー スペクトルを返します。

• x がベクトルまたはデータのベクトルをもつ timetable の場合、単一チャネルとして扱われます。

• x が 1 つの行列、行列変数をもつ timetable、または複数のベクトル変数をもつ timetable の場合、スペクトルは各チャネルに対して個別に計算され、p の個別の列に保存されます。

p = pspectrum(x,fs) は、fs のレートでサンプリングされたベクトルまたは行列信号のパワー スペクトルを返します。

p = pspectrum(x,t) は、t で指定された時点でサンプリングされたベクトルまたは行列信号のパワー スペクトルを返します。

p = pspectrum(___,type) は、関数で実行されるスペクトル解析の種類を指定します。type は、'power'、'spectrogram'、または 'persistence' のいずれかに指定します。この構文には、前の構文の入力引数を任意に組み合わせて含めることができます。

p = pspectrum(___,Name,Value) は、名前と値のペアの引数を使用して追加オプションを指定します。オプションには、周波数分解能帯域幅や、隣り合ったセグメント間のオーバーラップ率があります。

[p,f] = pspectrum(___) はスペクトル推定に対応する周波数を p に含めて返します。

[p,f,t] = pspectrum(___,'spectrogram') は、短時間のパワー スペクトルの推定の計算に使用されるウィンドウ セグメントの中心に対応する時点のベクトルも返します。

[p,f,pwr] = pspectrum(___,'persistence') は、パーシステンス スペクトルに含まれる推定に対応するパワー値のベクトルも返します。

出力引数を設定せずに pspectrum(___) を使用すると、現在の Figure ウィンドウにスペクトル推定がプロットされます。プロットでは、この関数が 10 log₁₀(p) を使用して p を dB に変換します。

例

正弦波のパワー スペクトル

2 チャネルの複素正弦波の 128 サンプルを生成します。

• 最初のチャネルは、単位振幅および $\pi /4$ ラジアン/サンプルの正規化された正弦波周波数をもちます。

• 2 番目チャネルは、$1/\sqrt{2}$ の振幅および $\pi /2$ ラジアン/サンプルの正規化された周波数をもちます。

各チャネルのパワー スペクトルを計算してプロットします。$0.15\pi$ ラジアン/サンプルから $0.6\pi$ ラジアン/サンプルの周波数範囲を拡大表示します。pspectrum は、信号の周波数成分がビンに正確に収まる場合、そのビン内の振幅が信号の真の平均パワーになるようにスペクトルをスケールします。複素指数の場合、平均パワーは振幅の 2 乗です。信号の離散フーリエ変換を計算して確認します。詳細については、確定的周期信号のパワーの測定を参照してください。

N = 128;
x = [1 1/sqrt(2)].*exp(1j*pi./[4;2]*(0:N-1)).';

[p,f] = pspectrum(x);

plot(f/pi,p)
hold on
stem(0:2/N:2-1/N,abs(fft(x)/N).^2)
hold off
axis([0.15 0.6 0 1.1])
legend("Channel "+[1;2]+", "+["pspectrum" "fft"])
grid

1 kHz で 296 ms 間サンプリングされ、ホワイト ガウス ノイズに組み込まれた正弦波信号を生成します。200 Hz の正弦波周波数と 0.1² のノイズ分散を指定します。信号とその時間情報を MATLAB® timetable に保存します。

Fs = 1000;
t = (0:1/Fs:0.296)';
x = cos(2*pi*t*200)+0.1*randn(size(t));
xTable = timetable(seconds(t),x);

信号のパワー スペクトルを計算します。スペクトルをデシベル単位で表し、プロットします。

[pxx,f] = pspectrum(xTable);

plot(f,pow2db(pxx))
grid on
xlabel('Frequency (Hz)')
ylabel('Power Spectrum (dB)')
title('Default Frequency Resolution')

正弦波のパワー スペクトルを再計算しますが、今度は粗いほうの周波数分解能 25 Hz を使用します。出力引数なしで関数 pspectrum を使用してスペクトルをプロットします。

pspectrum(xTable,'FrequencyResolution',25)

両側スペクトル

3 kHz で 1 秒間サンプリングされた信号を生成します。信号は凸二次チャープで、測定中に周波数が 300 Hz から 1,300 Hz に増加します。チャープはホワイト ガウス ノイズに組み込まれます。

fs = 3000;
t = 0:1/fs:1-1/fs;

x1 = chirp(t,300,t(end),1300,'quadratic',0,'convex') + ...
    randn(size(t))/100;

箱型ウィンドウを使用して、信号の両側パワー スペクトルを計算してプロットします。実信号の場合、pspectrum は既定では片側スペクトルをプロットします。両側スペクトルをプロットするには、TwoSided を true に設定します。

pspectrum(x1,fs,'Leakage',1,'TwoSided',true)

持続時間とサンプル レートが同じ複素数値信号を生成します。信号は正弦関数的に変化する周波数成分をもつチャープであり、ホワイト ノイズに組み込まれています。信号のスペクトログラムを計算し、ウォーターフォール プロットとして表示します。複素数値信号の場合、スペクトログラムは既定では両側です。

x2 = exp(2j*pi*100*cos(2*pi*2*t)) + randn(size(t))/100;

[p,f,t] = pspectrum(x2,fs,'spectrogram');

waterfall(f,t,p')
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
wtf = gca;
wtf.XDir = 'reverse';
view([30 45])

ウィンドウの漏れとトーンの分解能

100Hz で 2 秒間サンプリングされた 2 チャネルの信号を生成します。

1. 最初のチャネルは、20Hz トーンと 21Hz トーンで構成されます。どちらのトーンにも単位振幅があります。

2. 2 番目のチャネルも 2 つのトーンをもちます。1 つのトーンは単位振幅および 20Hz の周波数をもちます。もう 1 つのトーンは 1/100 の振幅および 30Hz の周波数をもちます。

fs = 100;
t = (0:1/fs:2-1/fs)';

x = sin(2*pi*[20 20].*t) + [1 1/100].*sin(2*pi*[21 30].*t);

ホワイト ノイズに信号を組み込みます。40 dB の S/N 比を指定します。信号をプロットします。

x = x + randn(size(x)).*std(x)/db2mag(40);

plot(t,x)

2 つのチャネルのスペクトルを計算し、表示します。

pspectrum(x,t)

スペクトル漏れの既定値 0.5 は、分解能帯域幅の約 1.29 Hz に相当します。1 番目のチャネルの 2 つのトーンは、分解されていません。2 番目のチャネルの 30Hz トーンは、他のトーンよりもかなり弱いにもかかわらず表示されています。

0.85 に漏れを増やし、約 0.74 Hz の分解能と等価にします。2 番目のチャネルの弱いトーンは、はっきりと表示されます。

pspectrum(x,t,'Leakage',0.85)

漏れを最大値に増やします。分解能帯域幅はおよそ 0.5Hz です。1 番目のチャネルの 2 つのトーンは、分解されています。2 番目のチャネルの弱いトーンは、大きいウィンドウのサイドローブによってマスクされています。

pspectrum(x,t,'Leakage',1)

関数 spectrogram と pspectrum の比較

電圧制御発振器と 3 つのガウス原子で構成される信号を生成します。信号は ${\mathit{f}}_{\mathit{s}}=2$ kHz で 2 秒間サンプリングされます。

fs = 2000;
tx = 0:1/fs:2;
gaussFun = @(A,x,mu,f) exp(-(x-mu).^2/(2*0.03^2)).*sin(2*pi*f.*x)*A';
s = gaussFun([1 1 1],tx',[0.1 0.65 1],[2 6 2]*100)*1.5;
x = vco(chirp(tx+.1,0,tx(end),3).*exp(-2*(tx-1).^2),[0.1 0.4]*fs,fs);
x = s+x';

短時間フーリエ変換

関数 pspectrum を使用して STFT を計算します。

• ${\mathit{N}}_{\mathit{x}}$ サンプル信号を、$80/2000=40$ ミリ秒の時間分解能に対応する長さ $\mathit{M}=80$ のサンプルのセグメントに分割します。

• $\mathit{L}=16$ 個のサンプル、または隣接するセグメント間の 20% のオーバーラップを指定します。

• カイザー ウィンドウを使用して各セグメントにウィンドウを適用し、漏れには $\ell =0.7$ を指定します。

M = 80;
L = 16;
lk = 0.7;

[S,F,T] = pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk);

関数 spectrogram で得られた結果と比較します。

• ウィンドウの長さとオーバーラップをサンプル単位で直接指定します。

• pspectrum は常にカイザー ウィンドウを $\mathit{g}\left(\mathit{n}\right)$ として使用します。漏れ $\ell$ とウィンドウの形状係数 $\beta$ は、$\beta =40\times \left(1-\ell \right)$ の関係にあります。

• pspectrum は、離散フーリエ変換の計算に常に ${\mathit{N}}_{\mathrm{DFT}}=1024$ 個の点を使用します。両側または中央揃えの周波数範囲の変換を計算したい場合にこの数値を指定することができます。一方、実信号の既定値である片側変換の場合、spectrogram は $1024/2+1=513$ 個の点を使用します。あるいは、この例にあるように、変換を計算したい周波数のベクトルを指定することができます。

• 信号を $\mathit{k}=\lfloor \frac{{\mathit{N}}_{\mathit{x}}-\mathit{L}}{\mathit{M}-\mathit{L}}\rfloor$ 個のセグメントに厳密に分割できない場合、spectrogram は信号を切り捨て、pspectrum は信号をゼロでパディングして追加のセグメントを作成します。出力を等価にするには、最後のセグメントと時間ベクトルの最後の要素を削除します。

• spectrogram は、振幅の 2 乗をスペクトログラムとする STFT を返します。pspectrum は、あらかじめ係数 ${\sum }_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$ で除算してから 2 乗した各セグメントのパワー スペクトルを返します。

• 片側変換の場合、pspectrum はスペクトログラムに追加の係数 2 を追加します。

g = kaiser(M,40*(1-lk));

k = (length(x)-L)/(M-L);
if k~=floor(k)
    S = S(:,1:floor(k));
    T = T(1:floor(k));
end

[s,f,t] = spectrogram(x/sum(g)*sqrt(2),g,L,F,fs);

関数 waterplot を使用して、これら 2 つの関数によって計算されたスペクトログラムを表示します。

subplot(2,1,1)
waterplot(sqrt(S),F,T)
title("pspectrum")

subplot(2,1,2)
waterplot(s,f,t)
title("spectrogram")

maxd = max(max(abs(abs(s).^2-S)))

maxd = 
2.4419e-08

パワー スペクトルと簡易プロット

関数 spectrogram は、各セグメントのパワー スペクトルまたはパワー スペクトル密度に対応する 4 番目の引数をもちます。pspectrum の出力と同様、ps 引数はあらかじめ 2 乗されており、正規化係数 ${\sum }_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$ を含みます。実信号の片側スペクトログラムの場合も、追加の係数 2 を含める必要があります。関数のスケーリング引数を "power" に設定します。

[~,~,~,ps] = spectrogram(x*sqrt(2),g,L,F,fs,"power");

max(abs(S(:)-ps(:)))

ans = 
2.4419e-08

pspectrum と spectrogram は、出力引数なしで呼び出された場合に、信号のスペクトログラムをデシベル単位でプロットします。片側スペクトログラムの場合は係数 2 を含めます。両方のプロットで同じになるようにカラーマップを設定します。x の範囲を同じ値に設定して、pspectrum プロットの最後にある追加のセグメントが表示されるようにします。spectrogram プロットの y 軸に周波数を表示します。

subplot(2,1,1)
pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk)
title("pspectrum")
cc = clim;
xl = xlim;

subplot(2,1,2)
spectrogram(x*sqrt(2),g,L,F,fs,"power","yaxis")
title("spectrogram")
clim(cc)
xlim(xl)

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency (Hz)")
    ylabel("Time (s)")
end

過渡信号のパーシステンス スペクトル

広帯域信号に組み込まれた狭帯域の干渉信号を可視化します。

1 kHz で 500 秒間サンプリングされたチャープを生成します。チャープの周波数は、測定中に 180 Hz から 220 Hz に増加します。

fs = 1000;
t = (0:1/fs:500)';

x = chirp(t,180,t(end),220) + 0.15*randn(size(t));

この信号には 210 Hz の正弦波も含まれます。正弦波の振幅は 0.05 で、正弦波は信号の全持続時間の 1/6 の時間のみ存在します。

idx = floor(length(x)/6);
x(1:idx) = x(1:idx) + 0.05*cos(2*pi*t(1:idx)*210);

信号のスペクトログラムを計算します。周波数範囲を 100 Hz ～ 290 Hz に制限します。1 秒の時間分解能を指定します。両方の信号成分が表示されます。

pspectrum(x,fs,'spectrogram', ...
    'FrequencyLimits',[100 290],'TimeResolution',1)

信号のパワー スペクトルを計算します。弱い正弦波はチャープによって不明瞭になります。

pspectrum(x,fs,'FrequencyLimits',[100 290])

信号のパーシステンス スペクトルを計算します。今度は、両方の信号成分が明瞭に表示されます。

pspectrum(x,fs,'persistence', ...
    'FrequencyLimits',[100 290],'TimeResolution',1)

チャープのスペクトログラムと再代入したスペクトログラム

1 kHz で 2 秒間サンプリングされた 2 次チャープを生成します。チャープの初期周波数は 100 Hz で、t = 1 秒のとき 200 Hz まで増大します。関数 pspectrum の既定の設定を使用してスペクトログラムを計算します。関数 waterfall を使用して、スペクトログラムをプロットします。

fs = 1e3;
t = 0:1/fs:2;
y = chirp(t,100,1,200,"quadratic");

[sp,fp,tp] = pspectrum(y,fs,"spectrogram");

waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
ylabel("Time (s)")
xlabel("Frequency (Hz)")

再割り当てされたスペクトログラムを計算して表示します。

[sr,fr,tr] = pspectrum(y,fs,"spectrogram",Reassign=true);

waterfall(fr,tr,sr')
set(gca,XDir="reverse",View=[60 60])
ylabel("Time (s)")
xlabel("Frequency (Hz)")

0.2 秒の時間分解能を使用してスペクトログラムを再計算します。出力引数なしで関数 pspectrum を使用して結果を視覚化します。

pspectrum(y,fs,"spectrogram",TimeResolution=0.2)

同じ時間分解能を使用して再代入されたスペクトログラムを計算します。

pspectrum(y,fs,"spectrogram",TimeResolution=0.2,Reassign=true)

ダイヤル トーン信号のスペクトログラム

4 kHz でサンプリングされた信号を作成します。これはデジタル電話のすべてのキーを押下した場合と似ています。信号を MATLAB® の timetable として保存します。

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

ver = [697 770 852 941];
hor = [1209 1336 1477];

tones = [];

for k = 1:length(ver)
    for l = 1:length(hor)
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones = [tones;tone;zeros(size(tone))];
    end
end

% To hear, type soundsc(tones,fs)

S = timetable(seconds(0:length(tones)-1)'/fs,tones);

信号のスペクトログラムを計算します。0.5 秒の時間分解能を指定し、隣接するセグメント間のオーバーラップとしてゼロを入力します。0.85 の漏れを指定します。これは、ハン ウィンドウでデータをウィンドウ処理することとほぼ等価です。

pspectrum(S,'spectrogram', ...
    'TimeResolution',0.5,'OverlapPercent',0,'Leakage',0.85)

スペクトログラムに、各キーが 0.5 秒間押下され、キー間の無音の一時停止期間が 0.5 秒であることが示されます。最初のトーンの周波数成分は 697 Hz と 1209 Hz 付近に集中しており、DTMF 規格の数字 '1' に対応します。

入力引数

x — 入力信号

ベクトル | 行列 | timetable

入力信号。ベクトル、行列、または MATLAB^® timetable として指定します。

• x が timetable の場合、増加する有限の行時間を含んでいなければなりません。

  メモ

  timetable が欠損している場合や時間点が重複している場合、欠損または重複する時間および非等間隔の時間をもつ timetable の整理のヒントを使用して修正できます。

• x がマルチチャネル信号を表す timetable の場合は、行列を含む単一変数か、ベクトルで構成される複数の変数のどちらかをもたなければなりません。

x が不等間隔サンプルの場合、pspectrum は信号を等間隔グリッドに内挿して、スペクトル推定を計算します。関数は線形内挿を使用し、サンプル時間が隣接する時間点の差の中央値に等しいと仮定します。不等間隔サンプル信号がサポートされるには、時間間隔の中央値と時間間隔の平均値は以下に従わなければなりません。

$$\frac{1}{100}<\frac{\text{Median time interval}}{\text{Mean time interval}}<100.$$

例: cos(pi./[4;2]*(0:159))'+randn(160,2) は、ホワイト ノイズに組み込まれた正弦波で構成される 2 チャネル信号です。

例: timetable(seconds(0:4)',rand(5,2)) は 1 Hz で 4 秒間サンプリングされた 2 チャネルの確率変数を指定します。

例: timetable(seconds(0:4)',rand(5,1),rand(5,1)) は 1 Hz で 4 秒間サンプリングされた 2 チャネルの確率変数を指定します。

データ型: single | double
複素数のサポート: あり

fs — サンプル レート

2π (既定値) | 正の数値スカラー

サンプル レート。正の数値スカラーとして指定します。

t — 時間値

ベクトル | datetime 配列 | duration 配列 | duration スカラー

時間値。ベクトル、datetime、duration 配列、またはサンプル間の時間間隔を表す duration スカラーとして指定します。

例: seconds(0:1/100:1) は 100 Hz でサンプリングされた 1 秒間を表す duration 配列です。

例: seconds(1) は、連続する信号サンプル間の 1 秒間の時間差を表す duration スカラーです。

type — 計算するスペクトルのタイプ

'power' (既定値) | 'spectrogram' | 'persistence'

計算するスペクトルのタイプ。'power'、'spectrogram' または 'persistence' として指定します。

• 'power' — 入力のパワー スペクトルを計算します。定常信号の周波数成分を解析する場合に、このオプションを使用します。詳細については、スペクトルの計算を参照してください。

• 'spectrogram' — 入力のスペクトログラムを計算します。信号の周波数成分が時間の経過と共に変化する様子を解析する場合に、このオプションを使用します。詳細については、スペクトログラムの計算を参照してください。

• 'persistence' — 入力のパーシステンス パワー スペクトルを計算します。特定の周波数成分が信号に存在する瞬間を可視化する場合に、このオプションを使用します。詳細については、パーシステンス スペクトルの計算を参照してください。

メモ

'spectrogram' および 'persistence' オプションはマルチチャネル入力をサポートしません。

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、コンマを使用して名前と値をそれぞれ区切り、Name を引用符で囲みます。

例: 'Leakage',1,'Reassign',true,'MinThreshold',-35 は、箱型ウィンドウを使用してデータをウィンドウ処理し、再代入されたスペクトル推定を計算して、-35 dB より小さいすべての値をゼロに設定します。

FrequencyLimits — 周波数帯域の範囲

[0 fs/2] (既定値) | 2 要素の数値ベクトル

周波数帯域の範囲。'FrequencyLimits' と 2 要素数値ベクトルから構成されるコンマ区切りのペアとして指定します。

• 入力に時間情報が含まれる場合、周波数帯域は Hz 単位で表されます。

• 入力に時間情報が含まれない場合、周波数帯域は正規化された単位のラジアン/サンプルで表されます。

既定の設定では、pspectrum はナイキスト範囲全体のスペクトルを計算します。

• 指定された周波数帯域にナイキスト範囲外の領域が含まれる場合、pspectrum はその周波数帯域を切り捨てます。

• 指定された周波数帯域が完全にナイキスト範囲外の場合、pspectrum はエラーをスローします。

ナイキスト範囲の詳細については、スペクトルの計算を参照してください。

x のサンプリングが等間隔でない場合、pspectrum は信号を等間隔グリッドに線形内挿して、隣接する時間点の差の中央値の逆数と等価の有効なサンプル レートを定義します。'FrequencyLimits' を有効なサンプル レートで表します。

例: [0.2*pi 0.7*pi] は、時間情報をもたない信号の 0.2π ラジアン/サンプル ～ 0.7π ラジアン/サンプルのスペクトルを計算します。

FrequencyResolution — 周波数分解能帯域幅

実数の数値スカラー

周波数分解能帯域幅。'FrequencyResolution' および実数の数値スカラーで構成されるコンマ区切りのペアとして指定します。入力に時間情報が含まれる場合は Hz 単位、含まれない場合は正規化された単位のラジアン/サンプルで表します。この引数は、'TimeResolution' と同時に指定することはできません。この引数の既定値は、入力データのサイズによって異なります。詳細については、スペクトログラムの計算を参照してください。

例: pi/100 は、π/100 ラジアン/サンプルの周波数分解能を使用して、時間情報をもたない信号のスペクトルを計算します。

Leakage — スペクトル漏れ

0.5 (既定値) | 0 ～ 1 の実数の数値スカラー

スペクトル漏れ。'Leakage' および 0 ～ 1 の実数の数値スカラーで構成されるコンマ区切りのペアとして指定します。'Leakage' はメインローブの幅に相対するカイザー ウィンドウのサイドローブの減衰を制御し、分解能の改善と漏れの低減との間でバランスを取ります。

• 漏れの値が大きいと、近接したトーンが分解されますが、近傍の弱いトーンがマスクされます。

• 漏れの値が小さいと、より大きいトーンの近傍で小さなトーンが見つかりますが、近い周波数が固まって不鮮明になります。

例: 'Leakage',0 は、漏れを少なくしてスペクトル分解能の低下を最小にします。

例: 'Leakage',0.85 は、ハン ウィンドウを使用したデータのウィンドウ処理を近似します。

例: 'Leakage',1 は、箱型ウィンドウを使用したデータのウィンドウ処理と等価で、漏れは最大になりますが、スペクトル分解能が改善されます。

MinThreshold — 非ゼロ値の下限

-Inf (既定値) | 実数スカラー

非ゼロ値の下限。'MinThreshold' と実数スカラーで構成されるコンマ区切りのペアで指定します。pspectrum は、type 引数の値に応じて異なる 'MinThreshold' を実装します。

• 'power' または 'spectrogram' — pspectrum は、10 log₁₀(p) ≤ 'MinThreshold' がゼロとなるように p の要素を設定します。'MinThreshold' をデシベル単位で指定します。

• 'persistence' — pspectrum は、'MinThreshold' より小さい p の要素をゼロに設定します。'MinThreshold' には 0 ～ 100% を指定します。

NumPowerBins — パーシステンス スペクトルのパワー ビンの数

256 (既定値) | 20 から 1024 の整数

パーシステンス スペクトルのパワー ビンの数。'NumPowerBins' と、20 ～ 1024 の整数で構成されるコンマ区切りのペアとして指定します。

OverlapPercent — 隣り合ったセグメント間のオーバーラップ

[0, 100) の範囲にある実数スカラー

スペクトログラムまたはパーシステンス スペクトルの隣り合ったセグメント間のオーバーラップ。'OverlapPercent' と [0, 100) の範囲の実数スカラーで構成されるコンマ区切りのペアとして指定します。この引数の既定値は、スペクトル ウィンドウによって異なります。詳細については、スペクトログラムの計算を参照してください。

Reassign — 再割り当てオプション

false (既定値) | true

再割り当てオプション。'Reassign' と logical で構成されるコンマ区切りのペアとして指定します。このオプションを true に設定した場合、pspectrum は、時間と周波数の再割り当てを実行することで、スペクトル推定の局所化が鮮明になります。再割り当て手法では、読み取りと解釈の容易なピリオドグラムとスペクトログラムが作成されます。この手法では、各スペクトル推定はビンの幾何学的中心ではなく、そのビンのエネルギー中心に再代入されます。この手法により、チャープとインパルスの厳密な局所化が行われます。

TimeResolution — スペクトログラムまたはパーシステンス スペクトルの時間分解能

実数スカラー

スペクトログラムまたはパーシステンス スペクトルの時間分解能。'TimeResolution' および実数スカラーで構成されるコンマ区切りのペアとして指定します。入力に時間情報が含まれる場合は秒単位、含まれない場合は整数の数のサンプルで表します。この引数は、スペクトログラムまたはパーシステンス スペクトル推定を形成する短時間パワー スペクトルの計算に使用するセグメントの持続時間を制御します。'TimeResolution' は、'FrequencyResolution' と同時に指定することはできません。この引数の既定値は、入力データのサイズと、指定された場合は周波数分解能によって異なります。詳細については、スペクトログラムの計算を参照してください。

TwoSided — 両側スペクトル推定

false | true

両側スペクトル推定。'TwoSided' と logical 値で構成されるコンマ区切りのペアとして指定します。

• このオプションが true の場合、この関数は[–π, π]に対する中央揃えの両側スペクトル推定を計算します。入力に時間情報が含まれる場合、推定は[–f_s/2, f_s/2]に対して計算されます。ここで、f_s は有効なサンプル レートです。

• このオプションが false の場合、この関数はナイキスト範囲[0, π]に対する片側スペクトル推定を計算します。入力に時間情報が含まれる場合、推定は[0, f_s/2]に対して計算されます。ここで、f_s は有効なサンプル レートです。合計したパワーを保存するために、この関数は 0 とナイキスト周波数以外のすべての周波数でパワーを 2 倍にします。このオプションは、実信号でのみ有効です。

指定しない場合、'TwoSided' の既定値は実数入力信号の場合は false、複素入力信号の場合は true です。

出力引数

p — スペクトル

ベクトル | 行列

スペクトル。ベクトルまたは行列として返されます。スペクトルのタイプとサイズは、type 引数の値によって異なります。

• 'power' — p は、x の各チャネルのパワー スペクトル推定を含みます。この場合、p のサイズは N_f × N_ch であり、ここで、N_f は f の長さ、N_ch は x のチャネル数です。pspectrum は、信号の周波数成分がビンに正確に収まる場合、そのビン内の振幅が信号の真の平均パワーになるようにスペクトルをスケールします。たとえば、正弦波の平均パワーは、正弦波の振幅の 2 乗の半分です。詳細については、確定的周期信号のパワーの測定を参照してください。

• 'spectrogram' — p は、x の短時間の、時間が局所化されたパワー スペクトルの推定を含みます。この場合、p は、サイズ N_f × N_t に含まれます。ここで N_f は f の長さ、N_t は t の長さです。

• 'persistence' — p は、パーセントで表され、与えられた時間と周波数位置で信号が与えられたパワー レベルの成分をもつ確率を含みます。この場合、p は、サイズ N_pwr × N_f に含まれます。ここで N_pwr は pwr の長さ、N_f は f の長さです。

f — スペクトル周波数

ベクトル

スペクトル周波数。ベクトルとして返されます。入力信号に時間情報が含まれる場合、f は Hz 単位で表される周波数を含みます。入力信号に時間情報が含まれない場合、周波数は正規化された単位のラジアン/サンプルで表されます。

t — スペクトログラムの時間値

ベクトル | datetime 配列 | duration 配列

スペクトログラムの時間値。秒単位の時間値のベクトルまたは duration 配列として返されます。入力に時間情報が含まれない場合、t はサンプル数を含みます。t は、短時間のパワー スペクトル推定の計算に使用されるデータ セグメントの中心に対応する時間値を含みます。

• pspectrum への入力が timetable である場合、t の形式は、入力 timetable の時間値と同じです。

• pspectrum への入力が、数値、duration または datetime 配列で指定された一連の時点でサンプリングされた数値ベクトルである場合、t のタイプと形式は、入力時間値と同じです。

• pspectrum への入力が連続するサンプル間の指定された時間差をもつ数値ベクトルである場合、t は duration 配列です。

pwr — パーシステンス スペクトルのパワー値

ベクトル

パーシステンス スペクトルのパワー値。ベクトルとして返されます。

詳細

スペクトルの計算

信号のスペクトルを計算するため、pspectrum は信号の全長で達成可能なスペクトル分解能と、膨大な FFT を計算することから生じるパフォーマンス上の制限との間の妥協点を見つけます。

• この関数は、可能な場合、カイザー ウィンドウを使用して信号全体の単一の修正ピリオドグラムを計算します。

• 合理的な時間で単一の修正ピリオドグラムを計算できない場合は、ウェルチ ピリオドグラムを計算します。信号をオーバーラップ セグメントに分割し、各セグメントにカイザー ウィンドウを適用して、セグメントのピリオドグラムを平均します。

スペクトル ウィンドウ処理

実際の信号はいずれも有限の時間に対してのみ測定可能です。この事実は、信号が周期的または無限長のいずれかであることを前提とする、フーリエ解析に無視できない影響をもたらします。異なる信号サンプルに異なる重みを割り当てる "スペクトル ウィンドウ処理" は、有限サイズの影響をシステマチックに処理します。

信号をウィンドウ処理する最も簡単な方法は、測定間隔外は等しくゼロであり、すべてのサンプルが等しく重要であると仮定することです。この "箱型ウィンドウ" には、両端にスペクトル リンギングをもたらす不連続のジャンプがあります。その他のスペクトル ウィンドウはすべて両端が細くなっており、信号のエッジに近づくにつれて、より小さい重みをサンプルに割り当てることで、この影響を小さくします。

ウィンドウ処理には常に、分解能の改善と漏れの低減という相反する目的の間で妥協することが伴います。

• "分解能" は、信号エネルギーが周波数空間にどのように分布しているかを正確に把握するための能力です。理想的な分解能のスペクトル アナライザーは、信号内に存在する 2 つの異なるトーン (基本的な正弦波) を、周波数がどれほど近くても識別できます。量的には、この能力はウィンドウの変換のメインローブの幅に関係します。

• "漏れ" は、有限の信号で、すべての周波数成分が周波数範囲全体にわたってエネルギー量を投影することです。スペクトル内の漏れの量は、隣接する強いトーンの存在下で、ノイズから弱いトーンを検出する能力によって測定できます。量的には、この能力はウィンドウの周波数変換のサイドローブ レベルに関係します。

• スペクトルを正規化することにより、その帯域幅内の純音は、完全に中央に位置する場合に正確な振幅をもちます。

分解能を高くすると、漏れも高くなります。その逆も同様です。範囲の一方の端で、箱型ウィンドウのメインローブは可能な限り狭くなり、サイドローブは最も高くなります。このウィンドウは、近接したトーンが同様のエネルギー量を持つ場合は分解できますが、そうでない場合は、より弱いトーンを見つけることに失敗します。もう一方の端で、高いサイドローブ抑制を持つウィンドウに広いメインローブがあり、ここでは近い周波数がかたまって不鮮明になります。

pspectrum はカイザー ウィンドウを使用してウィンドウ処理を実行します。カイザー ウィンドウの場合、メインローブによって取得された信号エネルギーの比率は、調整可能な "形状係数" β に最大に依存します。pspectrum が使用する形状係数の範囲は、箱型ウィンドウに対応する β = 0 から、広いメインローブが倍精度で表現可能なスペクトル エネルギーを基本的にすべて取得する β = 40 までです。β ≈ 6 の中間値はハン ウィンドウをかなり正確に近似します。β を制御するには、名前と値のペア 'Leakage' を使用します。'Leakage' に ℓ を設定した場合、ℓ と β は β = 40(1 – ℓ) と関連しています。詳細については、kaiser を参照してください。

時間領域における β = 5.7 の 51 点ハン ウィンドウおよび 51 点カイザー ウィンドウ	周波数領域における β = 5.7 の 51 点ハン ウィンドウおよび 51 点カイザー ウィンドウ

パラメーターとアルゴリズムの選択

信号のスペクトルを計算するため、pspectrum は最初に 2 つのトーンが分解可能なままどのくらい近接できるかを測定する "分解能帯域幅" を決定します。分解能帯域幅には、次の理論値があります。

$${\text{RBW}}_{\text{theory}}=\frac{\text{ENBW}}{t_{\max }-t_{\min }}.$$

• t_max – t_min ("レコード長") は、選択した信号領域の時間領域の持続時間です。

• ENBW はスペクトル ウィンドウの "等価ノイズ帯域幅" です。詳細については、enbw を参照してください。

  ENBW を制御するには、名前と値のペア 'Leakage' を使用します。この引数の最小値は、β = 40 のカイザー ウィンドウに対応します。最大値は、β = 0 のカイザー ウィンドウに対応します。

しかし実際には、pspectrum は分解能を低くする可能性があります。分解能を低くすることで、合理的な時間でスペクトルを計算して、限りのあるピクセル数で表示できるようにします。これらの実用的な理由から、pspectrum が使用できる最も低い分解能帯域幅は次になります。

$${\text{RBW}}_{\text{performance}}=\frac{f_{\text{span}}}{4096-1},$$

ここで、f_span は 'FrequencyLimits' を使用して指定された周波数帯域の幅です。'FrequencyLimits' を指定しない場合、pspectrum では、サンプル レートを f_span として使用します。RBW_performance は調整できません。

信号のスペクトルを計算するため、関数は 2 つの値の大きい方を選択します。値は "ターゲット分解能帯域幅" と呼ばれます。

$$\text{RBW}=\max ({\text{RBW}}_{\text{theory}},{\text{RBW}}_{\text{performance}}).$$

• 分解能帯域幅が RBW_theory の場合、pspectrum は信号全体に対して単一の "修正ピリオドグラム" を計算します。関数は、カイザー ウィンドウと、名前と値のペア 'Leakage' で制御される形状係数を使用します。詳細については、periodogram を参照してください。

• 分解能帯域幅が RBW_performance の場合、pspectrum は信号に対して "ウェルチ ピリオドグラム" を計算します。関数は以下を実行します。

  1. 信号をオーバーラップ セグメントに分割する。

  2. カイザー ウィンドウと指定した形状係数を使用して、各セグメントを個別にウィンドウ処理する。

  3. すべてのセグメントのピリオドグラムを平均化する。

  ウェルチの手続きは、オーバーラップ セクションによって指定された信号の異なる "実現" を平均化し、ウィンドウを使用して冗長データを除去することでスペクトル推定の分散を軽減するために設計されています。詳細については、pwelch を参照してください。

  • 各セグメント (または同様にウィンドウ) の長さは、次を使用して計算されます。

    $$\text{Segment length}=\frac{f_{\text{Nyquist}}\times \text{ENBW}}{\text{RBW}},$$

    ここで、f_Nyquist は "ナイキスト周波数" です。(エイリアシングがない場合、ナイキスト周波数は、隣接する時間点の差の中央値の逆数として定義される有効なサンプル レートの 1/2 になります。"ナイキスト範囲" は、実信号の場合は [0, f_Nyquist]、複素信号の場合は [–f_Nyquist, f_Nyquist] です)。

  • 1 歩の長さは、初期推定値

    $$\text{Stride length}=\frac{\text{Window length}}{2\times \text{ENBW}-1},$$

    を調整して、最初のウィンドウが正確に最初のセグメントの最初のサンプルで開始され、最後のウィンドウが正確に最後のセグメントの最後のサンプルで終わるようにすることで求められます。

スペクトログラムの計算

非定常信号の時間依存スペクトルを計算するために、pspectrum は信号をオーバーラップしたセグメントに分割し、各セグメントにカイザー ウィンドウを適用して、短時間フーリエ変換を計算し、変換を連結して行列にします。詳細については、Signal Processing Toolbox を使用したスペクトログラムの計算を参照してください。

非定常信号は周波数成分が時間と共に変化する信号です。非定常信号の "スペクトログラム" は、周波数成分の時間発展の推定です。非定常信号のスペクトログラムを作成するために、pspectrum は以下の手順に従います。

1. 信号を等しい長さのセグメントに分割します。セグメントは十分に短くなければならず、信号の周波数成分がセグメント内で感知されるほど変化してはいけません。セグメントはオーバーラップしている場合も、そうでない場合もあります。

2. 各セグメントにウィンドウを適用してスペクトルを計算し、"短時間フーリエ変換" を求めます。

3. セグメントのスペクトルを使用してスペクトログラムを作成します。

   • 出力引数を設定して呼び出された場合、スペクトルを連結して行列にします。

   • 出力引数を設定せずに呼び出された場合、各スペクトルのパワーをセグメントごとにデシベル単位で表示します。振幅を振幅依存のカラーマップをもつイメージとして並べて表します。

関数が計算できるのは、単一チャネル信号のスペクトログラムのみです。

信号のセグメント分割

スペクトログラムを作成するために、まず、信号をオーバーラップの可能性のあるセグメントに分割します。関数 pspectrum で、名前と値のペアの引数 'TimeResolution' および 'OverlapPercent' を使用して、セグメントの長さおよび隣接するセグメント間のオーバーラップの量を制御することができます。長さとオーバーラップを指定しない場合、関数は、信号の全体長に基づいて長さを選択します。オーバーラップ率は

$$\left(1-\frac{\text{1}}{2\times \text{ENBW}-1}\right)\times 100,$$

で求められ、ここで ENBW はスペクトル ウィンドウの "等価ノイズ帯域" です。詳細については、enbw およびスペクトルの計算を参照してください。

指定された時間分解能

• 信号に時間情報がない場合は、時間分解能 (セグメント長) をサンプル単位で指定します。時間分解能は 1 以上で信号の長さ以下の整数でなければなりません。

  信号に時間情報がある場合は、時間分解能を秒単位で指定します。この関数は結果をサンプル数に変換し、サンプル数以下で 1 以上の最も近い整数に丸めます。時間分解能は信号の持続時間以下でなければなりません。

• オーバーラップをセグメント長のパーセンテージで指定します。関数は結果をサンプル数に変換し、サンプル数以下の最も近い整数に丸めます。

既定の時間分解能

時間分解能を指定しない場合、pspectrum は信号全体の長さを使用してセグメントの長さを選択します。関数は時間分解能を ⌈N/d⌉ サンプル数として設定します。ここで、⌈⌉ 記号は天井関数を表し、N は信号の長さ、d は N に依存する除数です。

信号の長さ (N)	除数 (d)	セグメント長
2 サンプル – 63 サンプル	2	1 サンプル – 32 サンプル
64 サンプル – 255 サンプル	8	8 サンプル – 32 サンプル
256 サンプル – 2047 サンプル	8	32 サンプル – 256 サンプル
2048 サンプル – 4095 サンプル	16	128 サンプル – 256 サンプル
4096 サンプル – 8191 サンプル	32	128 サンプル – 256 サンプル
8192 サンプル – 16383 サンプル	64	128 サンプル – 256 サンプル
16384 サンプル – N サンプル	128	128 サンプル – ⌈N / 128⌉ サンプル

隣接するセグメント間のオーバーラップを指定することもできます。オーバーラップを指定すると、セグメント数が変化します。信号の端点を越えるセグメントにはゼロが付加されます。

7 サンプルの信号 [s0 s1 s2 s3 s4 s5 s6] を考えてみましょう。⌈7/2⌉ = ⌈3.5⌉ = 4 であるため、オーバーラップがない場合、関数は信号を長さ 4 の 2 つのセグメントに分割します。オーバーラップが増加すると、セグメント数が変化します。

オーバーラップするサンプル数	結果のセグメント数
0	s0 s1 s2 s3 s4 s5 s6 0
1	s0 s1 s2 s3 s3 s4 s5 s6
2	s0 s1 s2 s3 s2 s3 s4 s5 s4 s5 s6 0
3	s0 s1 s2 s3 s1 s2 s3 s4 s2 s3 s4 s5 s3 s4 s5 s6

pspectrum は、最後のセグメントが信号の端点を越えた場合は、信号にゼロを付加します。この関数は、セグメントの中心に対応する時点のベクトル t を返します。

セグメントのウィンドウ処理とスペクトルの計算

pspectrum は信号をオーバーラップしたセグメントに分割した後で、各セグメントにカイザー ウィンドウを適用します。ウィンドウの形状係数 β により、漏れを 'Leakage' の名前と値のペアを使用して調整できます。次に関数は、各セグメントのスペクトルを計算し、スペクトルを連結してスペクトログラム行列にします。セグメント スペクトルを計算するには、pspectrum はスペクトルの計算で説明している手順に従います。ただし、分解能帯域幅の下限が次である場合を除きます。

$${\text{RBW}}_{\text{performance}}=4\times \frac{f_{\text{span}}}{1024-1}.$$

スペクトルのパワーの表示

出力引数を設定しないで呼び出された場合、関数はカラー バーと共に既定の MATLAB カラーマップを使用して、短時間フーリエ変換のパワーをデシベル単位で表示します。カラー バーはスペクトログラムのパワー範囲全体を含みます。

パーシステンス スペクトルの計算

"パーシステンス スペクトル" は、与えられた周波数が信号内に存在する時間の割合を示す時間-周波数領域の表示です。パーシステンス スペクトルはパワー周波数空間のヒストグラムです。信号が変化する中で特定の周波数が信号内に存在する時間が長ければ長いほど、その時間の割合は大きくなるため、表示内の色が明るく ("熱く") なります。パーシステンス スペクトルを使用して、他の信号の中に隠れている信号を識別します。

パーシステンス スペクトルを計算するために、pspectrum は以下のステップを実行します。

1. 指定された漏れ、時間分解能、およびオーバーラップを使用してスペクトログラムを計算します。詳細については、スペクトログラムの計算を参照してください。

2. パワーと周波数の値を 2 次元ビンに分割します (パワー ビンの数を指定するには、名前と値のペア 'NumPowerBins' を使用します)。

3. 各時間値について、パワー スペクトルの対数の二変量ストグラムを再計算します。その時点に信号エネルギーが存在する各パワー周波数ビンについて、対応する行列要素を 1 ずつ増やします。すべての時間値のヒストグラムの総和を求めます。

4. パワーと周波数に対して累積したヒストグラムをプロットします。色は、正規化された割合で表したヒストグラム カウントの対数に比例します。ゼロ値を表現するために、考えられる最小の振幅の 1/2 を使用します。

パワー スペクトル

ヒストグラム

累積ヒストグラム

参照

[1] harris, fredric j. “On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform.” Proceedings of the IEEE^®. Vol. 66, January 1978, pp. 51–83.

[2] Welch, Peter D. "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms." IEEE Transactions on Audio and Electroacoustics. Vol. 15, June 1967, pp. 70–73.

拡張機能

C/C++ コード生成
MATLAB® Coder™ を使用して C および C++ コードを生成します。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

pspectrum 関数は GPU 配列入力をサポートしますが、次の使用上の注意および制限があります。

• パーシステンス スペクトルはサポートされていません。

• 再割り当てされたスペクトルまたはスペクトログラムはサポートされていません。

詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2017b で導入

R2023b: コード生成での timetable の使用

関数 pspectrum は、コード生成で MATLAB timetable をサポートします。

R2023a: [プロットの作成] ライブ エディター タスクを使用した関数出力の可視化

[プロットの作成] ライブ エディター タスクを使用して、pspectrum の出力を対話的に可視化することができるようになりました。さまざまなグラフ タイプを選択し、オプションのパラメーターを設定できます。タスクは、ライブ スクリプトの一部となるコードも自動的に生成します。

参考

アプリ

• 信号アナライザー

関数

• periodogram | pwelch | spectrogram | stft

トピック

• Signal Processing Toolbox を使用したスペクトログラムの計算
• 時間-周波数ギャラリー
• 時間-周波数解析の実践的基礎

Teaching Resources

• Digital Signal Processing: Signals and Filter Design (MathWorks Teaching Resources)