Main Content

stftmag2sig

STFT 振幅からの信号の再構成

R2020b 以降

説明

x = stftmag2sig(s,nfft) は、再構成された時間領域の実信号 x を返します。これは、Griffin-Lim アルゴリズムに基づき短時間フーリエ変換 (STFT) 振幅 s から推定されたものです。関数は、s が離散フーリエ変換 (DFT) 長 nfft を使用して計算されたと仮定します。

x = stftmag2sig(s,nfft,fs) は、s がレート fs でサンプリングされたと仮定して、再構成された信号を返します。

x = stftmag2sig(s,nfft,ts) は、s がサンプル時間 ts でサンプリングされたと仮定して、再構成された信号を返します。

x = stftmag2sig(___,Name=Value) は、名前と値の引数を使用して追加オプションを指定します。オプションには、FFT ウィンドウや、初期位相を指定するメソッドなどが含まれます。これらの引数を前の入力構文のいずれかに追加できます。たとえば、FrequencyRange="onesided",InitializePhaseMethod="random" は、ランダムな初期位相の片側 STFT から信号が再構成されることを指定します。

[x,t,info] = stftmag2sig(___) は、信号が再構成される時点や、再構成プロセスについての情報を含む構造体も返します。

すべて折りたたむ

正規化周波数が π/60 ラジアン/サンプル、DC 値が 1 の正弦波サンプル 512 個について考えます。信号の STFT を計算します。

n = 512;

x = cos(pi/60*(0:n-1)')+1;

S = stft(x);

STFT の振幅から正弦波を再構成します。元の信号と再構成後の信号をプロットします。

xr = stftmag2sig(abs(S),size(S,1));

plot(x)
hold on
plot(xr,"--",LineWidth=2)
hold off
legend("Original","Reconstructed")

計算を繰り返しますが、今回は信号をゼロでパディングして、エッジの影響を減らします。

xz = circshift([x; zeros(n,1)],n/2);

Sz = stft(xz);
xr = stftmag2sig(abs(Sz),size(Sz,1));

xz = xz(n/2+(1:n));
xr = xr(n/2+(1:n));

plot(xz)
hold on
plot(xr,"--",LineWidth=2)
hold off
legend("Original","Reconstructed")

計算を繰り返しますが、今回は x を、長さが 2 倍の信号のセグメントだと仮定して、エッジの影響を減らします。

xx = cos(pi/60*(-n/2:n/2+n-1)')+1;

Sx = stft(xx);
xr = stftmag2sig(abs(Sx),size(Sx,1));

xx = xx(n/2+(1:n));
xr = xr(n/2+(1:n));

plot(xx)
hold on
plot(xr,"--",LineWidth=2)
hold off
legend("Original","Reconstructed")

減少する 2 つのチャープと広帯域のスプラッター音を含むオーディオ信号を読み込みます。信号は 8192 Hz でサンプリングされています。信号の STFT をプロットします。波形を 128 サンプルのセグメントに分割し、ハミング ウィンドウを使用してセグメントにウィンドウを適用します。隣り合ったセグメント間のオーバーラップを 64 サンプル、DFT 点を 1024 に指定します。

load splat
ty = (0:length(y)-1)/Fs;

% To hear, type sound(y,Fs)

wind = hamming(128);
olen = 64;
nfft = 1024;

stft(y,Fs,Window=wind,OverlapLength=olen,FFTLength=nfft)

STFT の振幅と位相を計算します。

s = stft(y,Fs,Window=wind,OverlapLength=olen,FFTLength=nfft);

smag = abs(s);
sphs = angle(s);

STFT の振幅に基づいて信号を再構成します。STFT を計算したときと同じパラメーターを使用します。既定では、stftmag2sig は、位相をゼロに初期化し、100 回の最適化反復を使用します。

[x,tx,info] = stftmag2sig(smag,nfft,Fs, ...
    Window=wind,OverlapLength=olen);

% To hear the reconstruction, type sound(x,Fs)

元の信号と再構成後の信号をプロットします。より良く比較するために、再構成された信号を垂直方向と右方向にシフトします。

plot(ty,y,tx+500/Fs,x+1)
legend("Original","Reconstructed",Location="best")

最後の 2 つの繰り返し間の収束に向けた相対的改善を出力します。

impr = info.Inconsistency
impr = 0.0579

最適化反復の数を 2 倍にし、初期位相を STFT から実際の位相に設定して、再構成を改善します。元の信号と再構成後の信号をプロットします。より良く比較するために、再構成された信号の負の値をプロットし、それを垂直方向と右方向にシフトします。

[x,tx,info] = stftmag2sig(smag,nfft,Fs, ...
    Window=wind,OverlapLength=olen, ...
    MaxIterations=200,InitialPhase=sphs);

% To hear the reconstruction, type sound(x,Fs)

plot(ty,y,tx+500/Fs,-x+1)
legend("Original","Reconstructed",Location="best")

最後の 2 つの繰り返し間の収束に向けた相対的改善を出力します。

impr = info.Inconsistency
impr = 2.0848e-16

5 kHz で 4 秒間サンプリングされた信号を生成します。信号は、振動する振幅の領域および増加トレンドで変動する周波数の領域によって分離された、持続時間が減少していく一連のパルスで構成されています。

fs = 5000;
t = 0:1/fs:4-1/fs;

x = 10*besselj(0,1000*(sin(2*pi*(t+2).^3/60).^5));

信号の短時間フーリエ変換を計算してプロットします。信号を 128 サンプルのセグメントに分割し、ハン ウィンドウを使用してセグメントをウィンドウ処理します。隣り合ったセグメント間のオーバーラップを 96 サンプル、DFT 点を 128 に指定します。

win = hann(128);
olen = 96;
nfft = 128;
stft(x,fs,Window=win,OverlapLength=olen,FFTLength=nfft)

STFT の振幅と位相を計算します。

s = stft(x,fs,Window=win,OverlapLength=olen,FFTLength=nfft);

smag = abs(s);
sphs = angle(s);

STFT の振幅に基づいて信号を再構成します。STFT を計算したときと同じパラメーターを使用します。既定の Griffin-Lim アルゴリズムを使用する代わりに、ADAM オプティマイザーで勾配降下法を使用します。

[xrec,trec,info] = stftmag2sig(smag,nfft,fs, ...
    Window=win,Method="gd",Optimizer="adam");

元の信号と再構成後の信号をプロットします。より良く比較するために、再構成された信号を垂直方向にシフトして、両方が表示されるようにします。

plot(t,x,trec,xrec+12)
legend("Original","Reconstructed",Location="northwest")
ylim([-5 27])

元の信号と再構成信号の間の平均二乗誤差を計算します。

sum((x-xrec').^2)/length(x)
ans = 1.3620

再現性を得るために、乱数シードを既定値に設定します。2 チャネル正弦波の 160 サンプルを生成します。各チャネルにノイズを付加します。

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

  • 2 番目のチャネルは、単位振幅および π/2 ラジアン/サンプルの正規化された正弦波周波数をもちます。

信号の STFT を計算します。信号を 32 サンプルのセグメントに分割し、ハン ウィンドウを使用してセグメントをウィンドウ処理します。

rng default
N = 160;
t = 0:N-1;
win = hann(32);
signal = cos(pi./[4;2]*t)'+randn(N,2)/5;

s = stft(signal,Window=win);

STFT の振幅から信号を再構成します。L-BFGS オプティマイザーで勾配降下法を使用します。L-BFGS オプティマイザーは、各反復のステップ サイズを自動的に選択します。最適化プロセス中の各チャネルの正規化された不整合を追跡します。

[x,tx,info] = stftmag2sig(abs(s), ...
    size(s,1),Window=win, ...
    Method="gd",Optimizer="lbfgs", ...
    Display=true);
#Iteration  |  Normalized Inconsistency  
      1     |  1.2078e+00 1.2042e+00 
     20     |  3.2348e-02 9.6221e-03 
     40     |  1.9969e-02 6.8350e-03 
     60     |  2.5463e-03 1.2099e-03 
     80     |  2.4790e-03 8.4016e-04 
    100     |  5.4774e-03 7.1804e-04 

Decomposition process stopped.
The number of iterations reached the specified "MaxIterations" value of 100.

元の信号と再構成後の信号をプロットします。

tiledlayout(2,1)
nexttile
plot(t,signal(:,1))
hold on
plot(tx,x(:,1),"--",LineWidth=2)
hold off
title("First Channel")
legend("Original","Reconstructed", ...
    Location="northoutside")
nexttile
plot(t,signal(:,2))
hold on
plot(tx,x(:,2),"--",LineWidth=2)
hold off
title("Second Channel")

最後の 2 つの繰り返し間の収束に向けた相対的改善を出力します。

info.Inconsistency
ans = 1×2

    0.0055    0.0007

入力引数

すべて折りたたむ

STFT 振幅。行列または 3 次元配列として指定します。s は実数値の信号に対応しなければなりません。

  • s が単一チャネル信号に対応する場合、s を行列として指定します。

  • s がマルチチャネル信号に対応する場合、s を 3 次元配列として指定します。ここで、3 番目の次元はチャネルに対応します。

例: abs(stft(sin(pi/2*(0:255)),FFTLength=128)) は、正弦波の STFT 振幅を指定します。

例: abs(stft(chirp(0:1/1e3:1,25,1,50))) は、1 kHz でサンプリングされたチャープの STFT 振幅を指定します。

例: abs(stft(cos(pi./[4;2]*(0:159))'+randn(160,2)/5)) は、ノイズを含む 2 チャネル正弦波の STFT 振幅を指定します。

データ型: single | double

DFT 点の数。正の整数スカラーとして指定します。この引数が常に必要です。

データ型: single | double

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

サンプル時間。duration スカラーで指定します。ts の指定は、サンプル レート fs = 1/ts の設定と等価です。

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

データ型: duration

名前と値の引数

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

例: FrequencyRange="onesided",InitializePhaseMethod="random",Method="gd" は、勾配降下法を使用し、ランダムな初期位相の片側 STFT から信号が再構成されることを指定します。

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

例: "FrequencyRange","onesided","InitializePhaseMethod","random" は、ランダムな初期位相の片側 STFT から信号が再構成されることを指定します。

不整合の表示オプション。数値または logical の 1 (true) または 0 (false) として指定します。このオプションを true に設定すると、stftmag2sig は 20 回の最適化反復ごとに、正規化された不整合を表示します。関数は実行の最後に停止情報も表示します。

データ型: logical

STFT 振幅の周波数範囲。次のいずれかとして指定します。

  • "centered"s を中央揃えの両側 STFT の振幅として扱います。nfft が偶数の場合、s は区間 (–π, π] ラジアン/サンプルで計算されています。nfft が奇数の場合、s は区間 (–π, π) ラジアン/サンプルで計算されています。時間情報を指定すると、計算区間はそれぞれ (–fs, fs/2] サイクル/単位時間、(–fs, fs/2) サイクル/単位時間となります。ここで、fs はサンプル レートです。

  • "twosided"s を区間 [0, 2π) ラジアン/サンプルで計算された両側 STFT の振幅として扱います。時間情報を指定した場合、計算区間は [0, fs) サイクル/単位時間となります。

  • "onesided"s を片側 STFT の振幅として扱います。nfft が偶数の場合、s は区間 [0, π] ラジアン/サンプルで計算されています。nfft が奇数の場合、s は区間 [0, π) ラジアン/サンプルで計算されています。時間情報を指定すると、計算区間はそれぞれ [0, fs/2] サイクル/単位時間、[0, fs/2) サイクル/単位時間となります。ここで、fs はサンプル レートです。

データ型: char | string

再構成プロセスの不整合許容誤差。正のスカラーとして指定します。正規化された不整合が許容誤差よりも低い場合、再構成プロセスは停止します。

データ型: single | double

位相の初期化。"zeros" または "random" として指定します。InitializePhaseMethodInitialPhase のいずれか 1 つのみを指定します。

  • "zeros" — 関数は、位相をゼロとして初期化します。

  • "random" — 関数は、区間 [–π, π] で一様分布した乱数として位相を初期化します。

データ型: char | string

初期位相。実数の行列または 3 次元配列として指定します。InitialPhase の要素は [–π, π] の範囲内になければなりません。InitialPhases と同じサイズでなければなりません。InitializePhaseMethodInitialPhase のいずれか 1 つのみを指定します。

例: angle(stft(randn(1000,1))) は、ランダムな信号の短時間フーリエ変換の位相を指定します。

例: 2*pi*(rand(size(stft(randn(1000,1))))-1/2) は、区間 [–π, π] で一様分布したランダムな位相の行列を指定します。この行列のサイズは、ランダムな信号の短時間フーリエ変換と同じです。

.

データ型: single | double

入力時間次元。次のいずれかとして指定します。

  • "acrosscolumns" — 関数は s の時間次元が列に沿い、周波数次元が行に沿うと見なします。

  • "downrows" — 関数は s の時間次元が行に沿い、周波数次元が列に沿うと見なします。

データ型: char | string

最適化の最大反復回数。正の整数スカラーとして指定します。反復回数が MaxIterations を超えると、再構成プロセスは停止します。

データ型: single | double

信号の再構成アルゴリズム。次のいずれかとして指定します。

  • "gla" — Griffin と Lim によって提唱されたオリジナルの再構成アルゴリズム ([1]を参照)。

  • "fgla" — Perraudin、Balazs、Søndergaard によって提唱された高速 Griffin-Lim アルゴリズム ([2]を参照)。

  • "legla" — Le Roux、Kameoka、Ono、Sagayama によって提唱された高速アルゴリズム ([3]を参照)。

  • "gd" — Ji および Tie によって提唱された勾配降下法 ([4]を参照)。 (R2023b 以降)

メモ

勾配降下法を指定するには、Deep Learning Toolbox™ のライセンスを所有していなければなりません。

データ型: char | string

R2023b 以降

勾配降下法で使用するオプティマイザー。次のいずれかとして指定します。

  • "sgdm" — モーメンタム項付き確率的勾配降下法 (SGDM) オプティマイザー

  • "adam" — 適応モーメント推定 (ADAM) オプティマイザー

  • "rmsprop" —平方根平均二乗伝播 (RMSProp) オプティマイザー

  • "lbfgs" — メモリ制限 Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) オプティマイザー

この引数は、Method"gd" に設定されている場合にのみ適用されます。

データ型: char | string

隣り合ったセグメント間でオーバーラップするサンプルの数。Window の長さより小さい正の整数として指定します。信号を正常に再構成させるには、OverlapLength を、STFT 振幅の生成に使用されるオーバーラップしたセグメントの数と一致させる必要があります。OverlapLength は、省略するか空として指定すると、ウィンドウの長さの 75% 以下となる最大の整数に設定されます。これは、既定のハン ウィンドウの 96 サンプルです。

データ型: double | single

R2023b 以降

勾配降下法で使用される反復更新のステップ サイズ。正のスカラーとして指定します。この引数は、Method"gd" に設定され、Optimizer"lbfgs" ではない場合にのみ有効です。

データ型: single | double

"legla" の更新規則の打ち切り次数。正の整数として指定します。この引数は、Method"legla" に設定されている場合にのみ適用され、そのメソッドの各反復で更新される位相値の数を制御します。指定しない場合、TruncationOrder" は適応アルゴリズムを使用して決定されます。

データ型: single | double

高速 Griffin-Lim アルゴリズムの更新パラメーター。正のスカラーとして指定します。この引数は、Method"fgla" に設定されている場合にのみ適用され、そのメソッドの更新規則に対するパラメーターを指定します。

データ型: single | double

スペクトル ウィンドウ。ベクトルとして指定します。信号を正常に再構成させるには、Window を、STFT 振幅の生成に使用されるウィンドウと一致させる必要があります。ウィンドウを指定しない場合、またはウィンドウを空として指定する場合、関数は長さが 128 の周期的ハン ウィンドウを使用します。Window の長さは 2 以上でなければなりません。

利用可能なウィンドウのリストについては、ウィンドウを参照してください。

例: hann(128,"periodic")(1-cos(2*pi*(128:-1:1)'/128))/2 は両方とも、stftmag2sig で使用される既定のウィンドウを指定します。

データ型: double | single

出力引数

すべて折りたたむ

再構成された時間領域信号。ベクトルまたは行列として返されます。

メモ

xs を同じ長さになるようにする場合、(length(x)-noverlap)/(length(window)-noverlap) の値は整数でなければなりません。Window を使用して window の長さを指定し、OverlapLength を使用して noverlap を指定します。

信号が再構成される時点。ベクトルとして返されます。

再構成プロセス情報。次のフィールドを含む構造体として返されます。

  • ExitFlag — 終了フラグ。

    • 値が 0 であれば、最大反復回数に到達してアルゴリズムが停止したことを示します。

    • 値が 1 であれば、相対許容誤差を満たしてアルゴリズムが停止したことを示します。

  • NumIterations — 合計反復回数。

  • Inconsistency — 最後の 2 つの繰り返し間の収束に向けた相対的改善の平均。

  • ReconstructedPhase — 最後の反復における再構成された位相。

  • ReconstructedSTFT — 最後の反復における再構成された短時間フーリエ変換。

詳細

すべて折りたたむ

短時間フーリエ変換

短時間フーリエ変換 (STFT) を使用して、非定常信号の周波数成分が時間の経過と共に変化する様子を解析します。STFT の振幅の 2 乗は、信号の "スペクトログラム" 時間-周波数表現と呼ばれます。スペクトログラムの詳細や、関数 Signal Processing Toolbox™ を使用したスペクトログラムの計算方法については、Signal Processing Toolbox を使用したスペクトログラムの計算を参照してください。

信号の 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.

  • 短時間フーリエ変換は可逆変換です。逆変換プロセスは、ウィンドウが適用されたセグメントのオーバーラップ加算により、ウィンドウ エッジでの信号の減衰を補正します。詳細については、逆短時間フーリエ変換を参照してください。

  • 関数 istft は、信号の STFT の逆変換を行います。

  • 特定の条件下では、信号の "完全再構成" を実現することができます。詳細については、完全再構成を参照してください。

  • 関数 stftmag2sig は、STFT の振幅から再構成された信号の推定を返します。

正規化された不整合

"正規化された不整合" は、連続する最適化反復における再構成プロセスの収束に向けた改善を測定します。

正規化された不整合の定義は、信号の再構成アルゴリズム Method によって異なります。

  • Method"gla""fgla"、または "legla" の場合、不整合は次のように定義されます。

    Inconsistency=STFT(ISTFT(sest))sestsest,

    ここで、sest は各反復で推定された複素短時間フーリエ変換であり、双柱は行列ノルム、STFT は短時間フーリエ変換、ISTFT はその逆変換を表します。

  • Method"gd" の場合、不整合は次のように定義されます。

    Inconsistency=sest,isest,i1s,

    ここで、sest,i は i 回目の反復で推定された推定複素短時間フーリエ変換であり、s は入力短時間フーリエ変換の振幅です。

stftmag2sig は、MATLAB® 関数 norm を使用して行列ノルムを計算します。STFT とその逆の詳細については、短時間フーリエ変換逆短時間フーリエ変換を参照してください。

ヒント

  • 勾配降下法を使用していて再構成が十分でない場合、Displaytrue に設定します。反復中の不整合を観察します。不整合が減少しない場合は、再構成を改善するために StepSize を減らします。

  • 勾配降下法を使用している場合、通常は L-BFGS オプティマイザーが最良の結果を提供します。このオプティマイザーは、各反復のステップ サイズを自動的に選択します。ただし、L-BFGS オプティマイザーは、同じ回数の反復を実行するために、他のオプティマイザーよりも多くの計算時間を必要とする場合があります。

参照

[1] Griffin, Daniel W., and Jae S. Lim. "Signal Estimation from Modified Short-Time Fourier Transform." IEEE Transactions on Acoustics, Speech, and Signal Processing. Vol. 32, Number 2, April 1984, pp. 236–243. https://doi.org/10.1109/TASSP.1984.1164317.

[2] Perraudin, Nathanaël, Peter Balazs, and Peter L. Søndergaard. "A Fast Griffin-Lim Algorithm." In 2013 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, New Paltz, NY, October 20–23, 2013. https://doi.org/10.1109/WASPAA.2013.6701851.

[3] Le Roux, Jonathan, Hirokazu Kameoka, Nobutaka Ono, and Shigeki Sagayama. "Fast Signal Reconstruction from Magnitude STFT Spectrogram Based on Spectrogram Consistency." In Proceedings of the 13th International Conference on Digital Audio Effects (DAFx-10), Graz, Austria, September 6–10, 2010.

[4] Ji, Li, and Zhou Tie. “On Gradient Descent Algorithm for Generalized Phase Retrieval Problem.” In 2016 IEEE 13th International Conference on Signal Processing (ICSP), 320–25. Chengdu, China: IEEE, 2016. https://doi.org/10.1109/ICSP.2016.7877848.

拡張機能

バージョン履歴

R2020b で導入

すべて展開する