Main Content

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

Signal Processing Toolbox を使用したスペクトログラムの計算

Signal Processing Toolbox™ には、非定常信号のスペクトログラムを計算する 3 つの関数が用意されています。各関数には、異なる入力引数、既定値、および出力があります。最適な選択は、特定の用途によって異なります。

"非定常" 信号は、周波数成分が時間の経過とともに変化する信号です。"短時間フーリエ変換" (STFT) は、この周波数成分が信号の変化に伴って変化する様子を解析するために使用します。STFT の振幅の 2 乗は、信号の "スペクトログラム" 時間-周波数表現と呼ばれます。

スペクトログラムは、時間-周波数表現の候補の 1 つにすぎません。Signal Processing Toolbox および Wavelet Toolbox™ で利用可能な他の時間-周波数表現の概要については、時間-周波数ギャラリーを参照してください。関数 periodogram を使用した定常信号の処理については、FFT を使用したパワー スペクトル密度推定を参照してください。

スペクトログラム計算用の関数

Signal Processing Toolbox には、スペクトログラムの計算に使用できる次の関数があります。

  • spectrogram — 最大限の柔軟性を目指して設計されています。STFT と、セグメントごとのパワー スペクトル密度またはパワー スペクトルを計算します。再割り当てをサポートします。

  • stft — 可逆性と最大限の制御を目指して設計されています。STFT を計算します。dlstft および stftLayer によって使用されます。マルチチャネル入力をサポートします。これに対応する関数 istft は、逆 STFT を計算します。

  • pspectrum — 使いやすさを目指して設計されています。パワー スペクトルを計算します。信号アナライザーによって生成された解析スクリプトで使用されます。定常信号のスペクトルとパーシステンス スペクトルを計算できます。再割り当てをサポートします。

関数 spectrogram は、以下の説明で参考として使用されます。

カテゴリパラメーター関数
spectrogramstftpspectrum
入力SignalNx 個の要素をもつベクトル
  • Nx 個の要素をもつベクトル

  • Nx 行の行列

  • Nx 行の timetable

  • Nx 個の要素をもつベクトル

  • Nx 行の timetable

ウィンドウ (g(n))

2 番目の位置引数

(既定: ハミング ウィンドウ)

名前と値の引数 Window

(既定: 周期的ハン ウィンドウ)

カイザー ウィンドウのみ
ウィンドウの長さ (M)

サンプル数として指定

(既定: ⌊Nx/4.5⌋)

サンプル数として指定

(既定: 128)

名前と値の引数 TimeResolution
漏れ (ℓ)
  • ウィンドウにより異なる

  • カイザー ウィンドウを使用する場合は、形状係数 β を使用して調整する

  • ウィンドウにより異なる

  • カイザー ウィンドウを使用する場合は、形状係数 β を使用して調整する

カイザー ウィンドウの形状係数 β に関連付けられた名前と値の引数 Leakage: ℓ = 1 – β/40
オーバーラップ (L)

3 番目の位置引数として指定されるサンプル数

(既定: ウィンドウの長さの 50%)

名前と値の引数 OverlapLength で指定されるサンプル数

(既定: ウィンドウの長さの 75%)

名前と値の引数 OverlapPercent で指定されるセグメント長の割合

(既定: (112×ENBW1)×100, ここで、ENBW はウィンドウの等価ノイズ帯域幅)

DFT 点の数 (NDFT)

4 番目の位置引数

(既定: max(256,2⌈log2M⌉))

名前と値の引数 FFTLength

(既定: 128)

常に 1024
時間情報5 番目の位置引数として指定されるサンプル レート2 番目の位置引数として指定されるサンプル レートまたは時間ベクトル2 番目の位置引数として指定されるサンプル レートまたは時間ベクトル
 
関数呼び出し
fs = 100;
x = exp(2j*pi*20* ...
    (0:1/fs:2-1/fs)); 
M = 200; lk = 0.5;
g = kaiser(M,40*(1-lk)); 
L = 100; Ndft = 1024;
sps = abs( ...
    spectrogram(x,g,L,     ...
    Ndft,fs,"centered")    ...
    ).^2;
sts = abs( ...
    stft(x,fs,Window=g,    ...
    OverlapLength=L,       ...
    FFTLength=Ndft)        ...
    ).^2;
pss = pspectrum(x,fs,      ...
    "spectrogram",         ...
    TimeResolution=M/fs,   ...
    Leakage=lk,            ...
    OverlapPercent=L/M*100 ...
    )*sum(g)^2;
 
簡易プロット

fs = 6e2;
ts = 0:1/fs:2.05;
x = vco(sin(2*pi*ts).* ...
    exp(-ts),[0.1 0.4]*fs,fs);
M = 32; lk = 0.9;
g = kaiser(M,40*(1-lk));
L = 22; Ndft = 1024;

  • spectrogram の場合、"power","yaxis" を追加

  • stft の場合、FrequencyRange="onesided" を追加

spectrogram function convenience plot.

stft function convenience plot

pspectrum function convenience plot

 
出力周波数範囲
  • 引数 freqrange を使用して次のように制御される:

    • "onesided"NDFT が偶数値の場合、周波数範囲はゼロ周波数とナイキスト周波数 fs/2 でクローズとなります。

      NDFT が奇数値の場合、周波数範囲はゼロ周波数でクローズ、fs/2 でオープンとなります。

      (実数値の信号の既定値。)

    • "twosided"NDFT のすべての値について、周波数範囲はゼロ周波数でクローズ、fs でオープンとなります。

    • "centered"NDFT が偶数値の場合、周波数範囲は –fs/2 でオープン、fs/2 でクローズとなります。

      NDFT が奇数値の場合、周波数範囲は両端でオープンとなります。

      (複素数値の信号の既定値。)

  • STFT とスペクトログラムの計算を行う周波数をユーザーがベクトルで指定できる。

名前と値の引数 FrequencyRange を使用して次のように制御される:

  • "onesided"spectrogram と同じ。

  • "twosided"spectrogram と同じ。

  • "centered"spectrogram と同じ。

    (実数値の信号と複素数値の信号の既定値。)

logical の名前と値の引数 TwoSided を使用して次のように制御される:

  • false — ゼロ周波数と fs/2 で範囲がクローズとなります。

    (実数値の信号の既定値。)

  • true–fs/2fs/2 で範囲がクローズとなります。

    (複素数値の信号の既定値。)

時間間隔
  • 最後の完全なウィンドウ セグメントの後で打ち切られた信号。

  • セグメント中心における時間値。

  • 最後の完全なウィンドウ セグメントの後で打ち切られた信号。

  • セグメント中心における時間値。

  • 最後の完全なウィンドウ セグメントより後がゼロ パディングされた信号。

  • セグメント中心における時間値。

正規化
  • 最初の出力引数は STFT。その振幅の 2 乗がスペクトログラム。

  • 4 番目の出力引数は振幅の 2 乗。スペクトログラムを取得するには、ng(n))2 で乗算し、"power" オプションを指定する。

最初の出力引数は STFT。その振幅の 2 乗がスペクトログラム。
  • 最初の出力引数は振幅の 2 乗。

  • スペクトログラムを取得するには、ng(n))2 で乗算する。

PSD とパワー スペクトル
  • spectrogram の 4 番目の出力引数には、セグメント パワー スペクトル密度またはセグメント パワー スペクトルが格納される。

  • スペクトログラムはウィンドウ要素の和の 2 乗とパワー スペクトルを乗算したものに等しい。

  • spectrumtype の引数:

    • "psd"ENBW で乗算してパワー スペクトルを取得

      (既定)

    • "power" — パワー スペクトル

出力引数は STFT。
  • 出力引数はパワー スペクトル

  • PSD を取得するには、ENBW で除算する

 

STFT とスペクトログラムの定義

信号の STFT は、信号上の長さ M の "解析ウィンドウ" g(n) をスライドして、ウィンドウが適用されたデータの各セグメントの離散フーリエ変換 (DFT) を計算することによって算出されます。ウィンドウは、R サンプルの間隔で元の信号を飛び越えます。これは、隣り合ったセグメント間の L = M – R 個のサンプルのオーバーラップに相当します。ほとんどのウィンドウ関数は、スペクトル リンギングを回避するためにエッジで小さくなります。ウィンドウが適用された各セグメントの DFT は、時間と周波数の各点の振幅と位相を含む複素数値行列に対して追加されます。STFT 行列は次の列数をもちます。

k=NxLML

ここで、Nx は信号 x(n) の長さです。⌊⌋ 記号は床関数を表します。行列内の行数は、中央変換および両側変換の場合は DFT 点の数である NDFT と同じで、実数値信号の片側変換の場合は NDFT/2 に近い奇数に等しくなります。

STFT 行列 X(f)=[X1(f)X2(f)X3(f)Xk(f)] の m 番目の列には、時間 mR 付近を中心としたウィンドウが適用されたデータの DFT が含まれます。

Xm(f)=n=x(n)g(nmR)ej2πfn.

関数 spectrogram と STFT 定義の比較

600 Hz で 2 秒間サンプリングされた複素数値凸 2 次チャープで構成される信号を生成します。チャープの初期周波数は 250 Hz、最終周波数は 50 Hz です。

fs = 6e2;
ts = 0:1/fs:2;
x = chirp(ts,250,ts(end),50,"quadratic",0,"convex","complex");

関数 spectrogram

関数 spectrogram を使用して信号の STFT を計算します。

  • 信号を、それぞれ M=49 サンプルの長さのセグメントに分割します。

  • 隣接するセグメント間に L=11 個のサンプルのオーバーラップを指定します。

  • 最後の短いセグメントを破棄します。

  • 各セグメントにバートレット ウィンドウを適用します。

  • NDFT=1024 点で各セグメントの離散フーリエ変換を評価します。既定では、複素数値信号の場合、spectrogram は両側変換を計算します。

M = 49;
L = 11;
g = bartlett(M);
Ndft = 1024;

[s,f,t] = spectrogram(x,g,L,Ndft,fs);

関数 waterplot を使用して、STFT の振幅の 2 乗として定義されるスペクトログラムを計算し、表示します。

waterplot(s,f,t)

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

STFT の定義

定義を使用して Nx サンプルの信号の STFT を計算します。信号を Nx-LM-L 個のオーバーラップ セグメントに分割します。各セグメントにウィンドウを適用し、それぞれについて、NDFT 点での離散フーリエ変換を評価します。

[segs,~] = buffer(1:length(x),M,L,"nodelay");

X = fft(x(segs).*g,Ndft);

STFT の時間範囲と周波数範囲を計算します。

  • 時間値を求めるには、時間ベクトルをオーバーラップ セグメントに分割します。時間値は、下端でオープンの区間として扱われる各セグメントの中間点です。

  • 周波数値を求めるには、ゼロ周波数でクローズ、下端でオープンとなるナイキスト区間を指定します。

tbuf = ts(segs);
tint = mean(tbuf(2:end,:));

fint = 0:fs/Ndft:fs-fs/Ndft;

spectrogram の出力を定義と比較します。スペクトログラムを表示します。

maxdiff = max(max(abs(s-X)))
maxdiff = 0
waterplot(X,fint,tint)

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

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

関数 spectrogramstft の比較

1.4 kHz で 2 秒間サンプリングされたチャープで構成される信号を生成します。チャープの周波数は、測定時間中に 600 Hz から 100 Hz に線形的に減少します。

fs = 1400;
x = chirp(0:1/fs:2,600,2,100);

stft 既定

関数 spectrogramstft を使用して信号の STFT を計算します。関数 stft の既定の値を使用します。

  • 信号を 128 サンプルのセグメントに分割し、各セグメントに周期的ハン ウィンドウを適用します。

  • 隣接するセグメント間に 96 個のサンプルのオーバーラップを指定します。この長さはウィンドウの長さの 75% と等価です。

  • 128 の DFT 点を指定し、STFT の中央を周波数ゼロに揃え、周波数は Hz 単位で表します。

2 つの結果が等価であることを確認します。

M = 128;
g = hann(M,"periodic");
L = 0.75*M;
Ndft = 128;

[sp,fp,tp] = spectrogram(x,g,L,Ndft,fs,"centered");

[s,f,t] = stft(x,fs);

dff = max(max(abs(sp-s)))
dff = 0

関数 mesh を使用して、2 つの出力をプロットします。

nexttile
mesh(tp,fp,abs(sp).^2)
title("spectrogram")
view(2), axis tight

nexttile
mesh(t,f,abs(s).^2)
title("stft")
view(2), axis tight

Figure contains 2 axes objects. Axes object 1 with title spectrogram contains an object of type surface. Axes object 2 with title stft contains an object of type surface.

spectrogram 既定

関数 spectrogram の既定の値を使用して計算を繰り返します。

  • 信号を M=Nx/4.5 の長さのセグメントに分割します。ここで、Nx は信号の長さです。各セグメントにハミング ウィンドウを適用します。

  • セグメント間で 50% のオーバーラップを指定します。

  • FFT を計算するには、max(256,2log2M) 点を使用します。正の正規化周波数のみ、スペクトログラムを計算します。

M = floor(length(x)/4.5);
g = hamming(M);
L = floor(M/2);
Ndft = max(256,2^nextpow2(M));

[sx,fx,tx] = spectrogram(x);

[st,ft,tt] = stft(x,Window=g,OverlapLength=L, ...
    FFTLength=Ndft,FrequencyRange="onesided");

dff = max(max(sx-st))
dff = 0

関数 waterplot を使用して、2 つの出力をプロットします。どちらの場合も、周波数軸を π で除算します。stft 出力の場合、サンプル数を有効なサンプル レート 2π で除算します。

figure
nexttile
waterplot(sx,fx/pi,tx)
title("spectrogram")

nexttile
waterplot(st,ft/pi,tt/(2*pi))
title("stft")

Figure contains 2 axes objects. Axes object 1 with title spectrogram, xlabel Frequency/\pi, ylabel Samples contains an object of type patch. Axes object 2 with title stft, xlabel Frequency/\pi, ylabel Samples contains an object of type patch.

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

関数 spectrogrampspectrum の比較

電圧制御発振器と 3 つのガウス原子で構成される信号を生成します。信号は fs=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 を計算します。

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

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

  • カイザー ウィンドウを使用して各セグメントにウィンドウを適用し、漏れには =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 は常にカイザー ウィンドウを g(n) として使用します。漏れ とウィンドウの形状係数 β は、β=40×(1-) の関係にあります。

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

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

  • spectrogram は、振幅の 2 乗をスペクトログラムとする STFT を返します。pspectrum は、あらかじめ係数 ng(n) で除算してから 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")

Figure contains 2 axes objects. Axes object 1 with title pspectrum, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title spectrogram, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

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

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

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

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

max(abs(S(:)-ps(:)))
ans = 2.4419e-08

pspectrumspectrogram は、出力引数なしで呼び出された場合に、信号のスペクトログラムをデシベル単位でプロットします。片側スペクトログラムの場合は係数 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)

Figure contains 2 axes objects. Axes object 1 with title pspectrum, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image. Axes object 2 with title spectrogram, xlabel Time (s), ylabel Frequency (kHz) contains an object of type image.

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

中央揃えのスペクトログラムと片側スペクトログラムの計算

2 kHz で 2 秒間サンプリングされた実数値のチャープで構成される信号を生成します。

fs = 2000;
tx = 0:1/fs:2;
x = vco(-chirp(tx,0,tx(end),2).*exp(-3*(tx-1).^2), ...
    [0.1 0.4]*fs,fs).*hann(length(tx))';

両側スペクトログラム

信号の両側 STFT を計算してプロットします。

  • 信号を、それぞれ M=73 サンプルの長さのセグメントに分割します。

  • 隣接するセグメント間に L=24 個のサンプルのオーバーラップを指定します。

  • 最後の短いセグメントを破棄します。

  • 各セグメントにフラットトップ ウィンドウを適用します。

  • NDFT=895 点で各セグメントの離散フーリエ変換を評価します。このとき、点の数が奇数であることに注意します。

M = 73;
L = 24;
g = flattopwin(M);
Ndft = 895;
neven = ~mod(Ndft,2);

[stwo,f,t] = spectrogram(x,g,L,Ndft,fs,"twosided");

出力引数なしで関数 spectrogram を使用し、両側スペクトログラムをプロットします。

spectrogram(x,g,L,Ndft,fs,"twosided","power","yaxis")

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Frequency (kHz) contains an object of type image.

定義を使用して両側スペクトログラムを計算します。隣接するセグメント間で L サンプルがオーバーラップする M サンプルのセグメントに信号を分割します。各セグメントにウィンドウを適用し、それぞれについて、NDFT 点での離散フーリエ変換を計算します。

[segs,~] = buffer(1:length(x),M,L,"nodelay");

Xtwo = fft(x(segs).*g,Ndft);

時間範囲と周波数範囲を計算します。

  • 時間値を求めるには、時間ベクトルをオーバーラップ セグメントに分割します。時間値は、下端でオープンの区間として扱われる各セグメントの中間点です。

  • 周波数値を求めるには、ゼロ周波数でクローズ、上端でオープンとなるナイキスト区間を指定します。

tbuf = tx(segs);
ttwo = mean(tbuf(2:end,:));

ftwo = 0:fs/Ndft:fs*(1-1/Ndft);

spectrogram の出力を定義と比較します。関数 waterplot を使用して、スペクトログラムを表示します。

diffs = [max(max(abs(stwo-Xtwo))) max(abs(f-ftwo')) max(abs(t-ttwo))]
diffs = 1×3
10-12 ×

         0    0.2274    0.0002

figure
nexttile
waterplot(Xtwo,ftwo,ttwo)
title("Two-Sided, Definition")

nexttile
waterplot(stwo,f,t)
title("Two-Sided, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title Two-Sided, Definition, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title Two-Sided, spectrogram Function, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

中央揃えのスペクトログラム

信号の中央揃えのスペクトログラムを計算します。

  • 両側 STFT で使用したのと同じ時間値を使用します。

  • 関数 fftshift を使用して、STFT のゼロ周波数成分をスペクトルの中央にシフトします。

  • NDFT が奇数の場合、周波数範囲は両側でオープンとなります。NDFT が偶数の場合、周波数範囲は下端でオープン、上端でクローズとなります。

出力を比較してスペクトログラムを表示します。

tcen = ttwo;

if ~neven
    Xcen = fftshift(Xtwo,1);
    fcen = -fs/2*(1-1/Ndft):fs/Ndft:fs/2;
else
    Xcen = fftshift(circshift(Xtwo,-1),1);
    fcen = (-fs/2*(1-1/Ndft):fs/Ndft:fs/2)+fs/Ndft/2;
end

[scen,f,t] = spectrogram(x,g,L,Ndft,fs,"centered");

diffs = [max(max(abs(scen-Xcen))) max(abs(f-fcen')) max(abs(t-tcen))]
diffs = 1×3
10-12 ×

         0    0.2274    0.0002

figure
nexttile
waterplot(Xcen,fcen,tcen)
title("Centered, Definition")

nexttile
waterplot(scen,f,t)
title("Centered, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title Centered, Definition, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title Centered, spectrogram Function, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

片側スペクトログラム

信号の片側スペクトログラムを計算します。

  • 両側 STFT で使用したのと同じ時間値を使用します。

  • NDFT が奇数の場合、片側 STFT は両側 STFT の最初の (NDFT+1)/2 行により構成されます。NDFT が偶数の場合、片側 STFT は両側 STFT の最初の NDFT/2+1 行により構成されます。

  • NDFT が奇数の場合、周波数範囲はゼロ周波数でクローズ、ナイキスト周波数でオープンとなります。NDFT が偶数の場合、周波数範囲は両側でクローズとなります。

出力を比較してスペクトログラムを表示します。実数値信号の場合、"onesided" 引数はオプションです。

tone = ttwo;

if ~neven
    Xone = Xtwo(1:(Ndft+1)/2,:);
else
    Xone = Xtwo(1:Ndft/2+1,:);
end

fone = 0:fs/Ndft:fs/2;

[sone,f,t] = spectrogram(x,g,L,Ndft,fs);

diffs = [max(max(abs(sone-Xone))) max(abs(f-fone')) max(abs(t-tone))]
diffs = 1×3
10-12 ×

         0    0.1137    0.0002

figure
nexttile
waterplot(Xone,fone,tone)
title("One-Sided, Definition")

nexttile
waterplot(sone,f,t)
title("One-Sided, spectrogram Function")

Figure contains 2 axes objects. Axes object 1 with title One-Sided, Definition, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch. Axes object 2 with title One-Sided, spectrogram Function, xlabel Frequency (Hz), ylabel Time (s) contains an object of type patch.

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

セグメントの PSD とパワー スペクトルの計算

関数 spectrogram は、4 番目の出力引数として、各セグメントのパワー スペクトル密度 (PSD) またはパワー スペクトルのいずれかを含む行列をもちます。パワー スペクトルは、PSD にウィンドウの等価ノイズ帯域幅 (ENBW) を乗算した値と等価です。

1 kHz で 1 秒間サンプリングされた対数チャープで構成される信号を生成します。チャープの初期周波数は 400 Hz で、測定の最後には 10 Hz に減少します。

fs = 1000;
tt = 0:1/fs:1-1/fs;
y = chirp(tt,400,tt(end),10,"logarithmic");

サンプル レートによる PSD とパワー スペクトルのセグメント化

信号を 102 サンプルのセグメントに分割し、各セグメントにハン ウィンドウを適用します。隣り合ったセグメント間のオーバーラップを 12 サンプル、DFT 点を 1024 に指定します。

M = 102;
g = hann(M);
L = 12;
Ndft = 1024;

既定の PSD スペクトル タイプで信号のスペクトログラムを計算します。STFT、およびセグメントのパワー スペクトル密度の配列を出力します。

[s,f,t,p] = spectrogram(y,g,L,Ndft,fs);

スペクトル タイプを "power" に指定して計算を繰り返します。STFT、およびセグメントのパワー スペクトルの配列を出力します。

[r,~,~,q] = spectrogram(y,g,L,Ndft,fs,"power");

どちらの場合もスペクトログラムが同じになることを確認します。周波数に対数スケールを使用して、スペクトログラムをプロットします。

max(max(abs(s).^2-abs(r).^2))
ans = 0
waterfall(f,t,abs(s)'.^2)
set(gca,XScale="log",...
    XDir="reverse",View=[30 50])

Figure contains an axes object. The axes object contains an object of type patch.

パワー スペクトルが、パワー スペクトル密度にウィンドウの ENBW を乗算した値と等価であることを確認します。

max(max(abs(q-p*enbw(g,fs))))
ans = 1.1102e-16

セグメントのパワー スペクトルの行列が、スペクトログラムに比例していることを確認します。比例係数は、ウィンドウ要素の和の 2 乗です。

max(max(abs(s).^2-q*sum(g)^2))
ans = 3.4694e-18

正規化周波数による PSD とパワー スペクトルのセグメント化

計算を繰り返しますが、今度は正規化周波数で計算します。サンプル レートを 2π として指定した場合、結果は同じになります。

[~,~,~,pn] = spectrogram(y,g,L,Ndft);
[~,~,~,qn] = spectrogram(y,g,L,Ndft,"power");

max(max(abs(qn-pn*enbw(g,2*pi))))
ans = 1.1102e-16

参考

アプリ

関数

関連するトピック

外部の Web サイト