istft

逆短時間フーリエ変換

構文

x = istft(s)

x = istft(s,fs)

x = istft(s,ts)

x = istft(___,Name,Value)

[x,t] = istft(___)

説明

x = istft(s) は s の 逆短時間フーリエ変換 (ISTFT) を返します。

x = istft(s,fs) は、サンプル レート fs を使用して s の ISTFT を返します。

x = istft(s,ts) は、サンプル時間 ts を使用して ISTFT を返します。

x = istft(___,Name,Value) は、名前と値のペアの引数を使用して追加オプションを指定します。オプションには、FFT ウィンドウの長さとオーバーラップされたサンプルの数が含まれます。これらの引数を前の入力構文のいずれかに追加できます。

[x,t] = istft(___) は、ISTFT が評価される信号時間を返します。

例

マルチチャネル信号の ISTFT

1 kHz で 1 秒間サンプリングされた 3 つの異なるチャープで構成される 3 チャネルの信号を生成します。

1. 最初のチャネルは凹二次チャープで構成され、瞬時周波数は t = 0 では 100 Hz であり、t = 1 秒では 300 Hz になります。初期位相は 45 度です。

2. 2 番目のチャネルは凸二次チャープで構成され、瞬時周波数は t = 0 では 200 Hz であり、t = 1 秒では 600 Hz になります。

3. 3 番目のチャネルは対数チャープで構成され、瞬時周波数は t = 0 では 300 Hz であり、t = 1 秒では 500 Hz になります。

長さが 256 の周期的ハミング ウィンドウと 15 サンプルのオーバーラップ長を使用して、マルチチャネル信号の STFT を計算します。

fs = 1e3;
t = 0:1/fs:1-1/fs;
x = [chirp(t,100,1,300,'quadratic',45,'concave');
      chirp(t,200,1,600,'quadratic',[],'convex');
      chirp(t,300,1,500,'logarithmic')]'; 
  
[S,F,T] = stft(x,fs,'Window',hamming(256,'periodic'),'OverlapLength',15);

最初と 2 番目のチャネルについて、元のバージョンと再構成されたバージョンをプロットします。

[ix,ti] = istft(S,fs,'Window',hamming(256,'periodic'),'OverlapLength',15);

plot(t,x(:,1)','LineWidth',1.5)
hold on
plot(ti,ix(:,1)','r--')
hold off
legend('Original Channel 1','Reconstructed Channel 1')

plot(t,x(:,2)','LineWidth',1.5)
hold on
plot(ti,ix(:,2)','r--')

legend('Original Channel 2','Reconstructed Channel 2')

異なる合成ウィンドウと解析ウィンドウをもつフェーズ ボコーダー

フェーズ ボコーダーは周波数領域でオーディオを変換することによってタイム ストレッチとピッチ スケーリングを実行します。この図は、フェーズ ボコーダーの実装に関連する操作を示します。

フェーズ ボコーダーは、ホップ サイズ ${\mathit{R}}_1$ の解析ウィンドウを持つ信号の STFT を実行してから、ホップ サイズ ${\mathit{R}}_2$ の合成ウィンドウを持つ ISTFT を実行します。ボコーダーは WOLA 方法を利用します。信号のタイム ストレッチを行うために、解析ウィンドウは合成よりも大きい数のオーバーラップ サンプルを使用します。結果として、周波数成分が同じだけ残っているにもかかわらず、入力時点よりも出力時点でのサンプルが多くなります (${\mathit{N}}_{\mathit{S},\mathrm{Out}}>{\mathit{N}}_{\mathit{S},\mathrm{In}}$)。高いサンプル レートで再生することでこの信号のピッチ スケーリングができます。つまり、元の持続時間のまま、より高いピッチで信号を作成します。

8192 Hz でサンプリングされたヘンデルの「ハレルヤ コーラス」の断章を含むオーディオ ファイルを読み込みます。

load handel

長さが 512 のルート-ハン ウィンドウを設計します。解析オーバーラップの長さを 192、合成オーバーラップの長さを 166 に設定します。

wlen = 512;
win = sqrt(hann(wlen,'periodic'));
noverlapA = 192;
noverlapS = 166;

オーバーラップ 192 の解析ウィンドウとオーバーラップ 166 の合成ウィンドウを使用して、フェーズ ボコーダーを実装します。

S = stft(y,Fs,'Window',win,'OverlapLength',noverlapA);
iy = istft(S,Fs,'Window',win,'OverlapLength',noverlapS);

%To hear, type soundsc(y,Fs), pause(10), soundsc(iy,Fs);

解析ウィンドウと合成ウィンドウが同じでオーバーラップの長さが変更される場合、調整が必要になるゲイン/損失が追加されます。これは、フェーズ ボコーダーを実装するための一般的な方法です。

ホップ率を計算して、それを使用して再構成後の信号のゲインを調整します。ホップ率を使用して、ピッチ シフトしたデータの周波数も計算します。

hopRatio = (wlen-noverlapS)/(wlen-noverlapA);
iyg = iy*hopRatio;
Fp = Fs*hopRatio;

%To hear, type soundsc(iyg,Fs), pause(15), soundsc(iyg,Fp);

元の信号と固定されたゲインでタイム ストレッチした信号をプロットします。

plot((0:length(iyg)-1)/Fs,iyg,(0:length(y)-1)/Fs,y)
xlabel('Time (s)')
xlim([0 (length(iyg)-1)/Fs])
legend('Time Stretched Signal with Fixed Gain','Original Signal','Location','best')

タイム ストレッチした信号とピッチ シフトした信号を同じプロットで比較します。

plot((0:length(iy)-1)/Fs,iy,(0:length(iy)-1)/Fp,iy)
xlabel('Time (s)')
xlim([0 (length(iyg)-1)/Fs])
legend('Time Stretched Signal','Pitch Shifted Signal','Location','best') 

ピッチ シフトするデータの効果をより理解するには、2 秒を超える次の周波数 Fs の正弦波について考えます。

t = 0:1/Fs:2;
x = sin(2*pi*10*t);

オーバーラップの長さが 192 の短時間フーリエ変換と 166 の逆短時間フーリエ変換を計算します。

Sx = stft(x,Fs,'Window',win,'OverlapLength',noverlapA);
ix = istft(Sx,Fs,'Window',win,'OverlapLength',noverlapS);

1 つのプロットに元の信号、別のプロットにタイム ストレッチした信号とピッチ シフトした信号をプロットします。

tiledlayout(2,1)
nexttile
plot((0:length(ix)-1)/Fs,ix,'LineWidth',2) 
xlabel('Time (s)')
ylabel('Signal Amplitude')
xlim([0 (length(ix)-1)/Fs])
legend('Time Stretched Signal') 
nexttile
plot((0:length(x)-1)/Fs,x)
hold on
plot((0:length(ix)-1)/Fp,ix,'--','LineWidth',2)
legend('Original Signal','Pitch Shifted Signal','Location','best')
hold off
xlabel('Time (s)')
ylabel('Signal Amplitude')
xlim([0 (length(ix)-1)/Fs])

ゼロ パディングする複素信号の ISTFT

周波数 1 kHz、持続時間 2 秒間の複素正弦波を生成します。

fs = 1e3;
ts = 0:1/fs:2-1/fs;

x = exp(2j*pi*100*cos(2*pi*2*ts));

長さが 100 の周期的ハン ウィンドウを設計し、オーバーラップ サンプルの数を 75 に設定します。COLA 準拠に対するウィンドウとオーバーラップの長さを確認します。

nwin = 100;
win = hann(nwin,'periodic');
noverlap = 75;

tf = iscola(win,noverlap)

tf = logical
   1

信号をゼロ パディングして、エッジの影響を除去します。切り捨てを避けるために、入力信号をゼロでパディングします。その結果、

$$\frac{length(xZero)-noverlap}{nwin-noverlap}$$

は整数になります。FFT 長を 128 に設定します。複素信号の短時間フーリエ変換を計算します。

nPad = 100;
xZero = [zeros(1,nPad) x zeros(1,nPad)];
fftlen = 128;
s = stft(xZero,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen);

逆短時間フーリエ変換を計算して、完全再構成のためにゼロを削除します。

[is,ti] = istft(s,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen);
is(1:nPad) = [];
is(end-nPad+1:end) = [];
ti = ti(1:end-2*nPad);

元の信号と再構成後の信号の実数部をプロットします。信号の虚数部も完全に再構成されます。

plot(ts,real(x))
hold on
plot(ti,real(is),'--')
xlim([0 0.5])
xlabel('Time (s)')
ylabel('Amplitude (V)')
legend('Original Signal','Reconstructed Signal') 
hold off

COLA 準拠のウィンドウとオーバーラップを使用した実信号の ISTFT

2 kHz で 1 秒間サンプリングされた正弦波を生成します。

fs = 2e3;
t = 0:1/fs:1-1/fs;
x = 5*sin(2*pi*10*t);

長さが 120 の周期的ハミング ウィンドウを設計します。オーバーラップが 80 サンプルのウィンドウの COLA 制約を確認します。ウィンドウとオーバーラップの組み合わせは COLA 準拠です。

win = hamming(120,'periodic');
noverlap = 80;
tf = iscola(win,noverlap)

tf = logical
   1

FFT 長を 512 に設定します。短時間フーリエ変換を計算します。

fftlen = 512;
s = stft(x,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen);

逆短時間フーリエ変換を計算します。

[X,T] = istft(s,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',fftlen,'Method','ola','ConjugateSymmetric',true);

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

plot(t,x,'b')
hold on
plot(T,X,'-.r')
xlabel('Time (s)')
ylabel('Amplitude (V)')
title('Original and Reconstructed Signal')
legend('Original Signal','Reconstructed Signal')
hold off

入力引数

s — 短時間フーリエ変換

行列 | 3 次元配列

短時間フーリエ変換。行列または 3 次元配列として指定します。単一チャネル信号の場合は、s を列方向に時間が増加し、行方向に周波数が増加する行列として指定します。マルチチャネル信号の場合は、s を 3 番目の次元がチャネルに対応する 3 次元配列として指定します。周波数ベクトルと時間ベクトルは stft の出力として取得されます。

メモ

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

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

fs — サンプル レート

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

サンプル レート (ヘルツ単位)。正のスカラーとして指定します。

データ型: double | single

ts — サンプル時間

duration スカラー

サンプル時間。duration スカラーで指定します。

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

データ型: duration

名前と値の引数

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

例: istft(s,Window=win,OverlapLength=50,FFTLength=128) は、隣接するセグメント間でのオーバーラップが 50 サンプル、DFT 点が 128 のウィンドウ win を使用してデータをウィンドウ処理します。

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

例: istft(s,'Window',win,'OverlapLength',50,'FFTLength',128) は、隣接するセグメント間でのオーバーラップが 50 サンプル、DFT 点が 128 のウィンドウ win を使用してデータをウィンドウ処理します。

Window — ウィンドウ処理関数

hann(128,"periodic") (既定値) | ベクトル

ウィンドウ処理関数。ベクトルとして指定します。ウィンドウを指定しない場合、またはウィンドウを空として指定する場合、関数は長さが 128 の周期的ハン ウィンドウを使用します。Window の長さは 2 以上でなければなりません。

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

例: hann(N+1) と (1-cos(2*pi*(0:N)'/N))/2 は、いずれも長さ N + 1 のハン ウィンドウを指定します。

データ型: double | single

OverlapLength — オーバーラップするサンプルの数

ウィンドウの長さの 75% (既定値) | 非負の整数

オーバーラップするサンプルの数。window の長さより小さな非負の整数として指定します。OverlapLength を省略するか、空として指定する場合、ウィンドウの長さの 75% より少ない最大整数に設定されます。これは、既定のハン ウィンドウの 96 サンプルになります。

データ型: double | single

FFTLength — DFT 点の数

128 (既定値) | 正の整数

DFT 点の数。ウィンドウの長さ以上の正の整数として指定します。完全な時間領域の再構成を行うには、stft で使用されるものに一致するように DFT 点の数を設定しなければなりません。

データ型: double | single

Method — オーバーラップ加算の方法

"wola" (既定値) | "ola"

オーバーラップ加算の方法。以下のいずれかとして指定します。

• "wola" — 重み付きオーバーラップ加算

• "ola" — オーバーラップ加算

ConjugateSymmetric — 元の信号の共役対称性

false (既定値) | true

元の信号の共役対称性。true または false として指定します。このオプションを true に設定する場合、istft は、入力 s が対称であると仮定します。そうでない場合は、対称性の仮定は作成されません。s が丸め誤差のために厳密に共役対称性でない場合、名前と値のペアの true への設定では STFT は共役対称として扱われます。s が共役対称である場合、逆変換の計算が速くなり、出力は実数になります。

FrequencyRange — STFT 周波数範囲

"centered" (既定値) | "twosided" | "onesided"

STFT 周波数範囲。"centered"、"twosided"、または "onesided" として指定します。

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

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

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

  メモ

  この引数が "onesided" に設定されている場合、istft は、正のナイキスト範囲の値が合計パワーを保存せずに計算されていると仮定します。

例については、STFT 周波数範囲を参照してください。

データ型: char | string

InputTimeDimension — 入力時間次元

"acrosscolumns" (既定値) | "downrows"

入力時間次元。"acrosscolumns" または "downrows" として指定します。この値が "downrows" に設定されている場合、istft は s の時間次元が行に沿い、周波数が列に沿うと見なします。この値が "acrosscolumns" に設定されている場合、関数 istft は s の時間次元が列に沿い、周波数次元が行に沿うと見なします。

出力引数

x — 再構成後の信号

ベクトル | 行列

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

データ型: single | double

t — 時点

ベクトル

時点。ベクトルとして返されます。

• サンプル レート fs が指定された場合、t には時間値 (秒単位) が含まれます。

• 持続時間 ts が指定された場合、t は、時間形式が入力持続時間と同じになり、duration 配列になります。

• 時間情報が指定されない場合、t はサンプル数を含みます。

データ型: double | single

詳細

逆短時間フーリエ変換

逆短時間フーリエ変換は、STFT の各 DFT ベクトルの IFFT を取り、逆変換した信号をオーバーラップ加算することで計算されます。ISTFT は次のように計算されます。

$$\begin{array}{c}x(n)={\displaystyle \underset{-1/2}{\overset{1/2}{\int }}{\displaystyle \sum _{m=-\infty }^{\infty }{X}_m(f)e^{j2\pi fn}}}df\\ ={\displaystyle \sum _{m=-\infty }^{\infty }{\displaystyle \underset{-1/2}{\overset{1/2}{\int }}{X}_m(f)e^{j2\pi fn}}}df\\ ={\displaystyle \sum _{m=-\infty }^{\infty }{x}_m}(n),\end{array}$$

ここで、$$R$$ は連続する DFT 間のホップ サイズで、$$X_m$$ は時間 $$mR,$$ と $$x_m(n)=x(n)g(n-mR)$$ 付近を中心としたウィンドウが適用されたデータの DFT です。逆 STFT は、"解析ウィンドウ" $$g(n)$$ が元の信号のウィンドウ処理に使用されて $$c$$ が定数の $${\displaystyle \sum _{m=-\infty }^{\infty }{g}^{a+1}(n-mR)=c}\forall n\in \mathbb{Z},$$ と同じ長さの元の信号の完全再構成です。次の図は、元の信号の再構成における各ステップを示したものです。

定数オーバーラップ加算 (COLA) 制約

変更していないスペクトルを正しく再構成させるには、解析ウィンドウが COLA 制約を満たさなければなりません。通常、解析ウィンドウが条件 $${\displaystyle \sum _{m=-\infty }^{\infty }{g}^{a+1}(n-mR)=c}\forall n\in \mathbb{Z},$$ を満たす場合、ウィンドウは COLA 準拠であると見なされます。さらに、COLA 準拠は弱いまたは強いのいずれかで説明されます。

• 弱い COLA 準拠とは、

  $$G(f_k)=0,k=1,2,\dots ,R-1,f_k\triangleq \frac{k}{R}.$$

  となるフレームレートの高調波で解析ウィンドウのフーリエ変換がゼロであることを意味します。

  エイリアスのキャンセルはスペクトル変更によって妨害されます。弱い COLA は、周波数領域のエイリアスのキャンセルに依存しています。したがって、信号に何らかのスペクトル変更が行われていなければ、弱い COLA 準拠ウィンドウを使用して完全な再構成が可能です。

• 強い COLA 準拠の場合、ウィンドウのフーリエ変換は

  $$G(f)=0,f\ge \frac{1}{2R}.$$

  となるフレーム レートによってダウンサンプリングしながら常に帯域制限しなければなりません。この方程式は、強い COLA 制約によって許可されるエイリアシングがないことを示します。さらに、強い COLA 準拠の場合、定数 $$c$$ の値は、1 と等しくなければなりません。通常、短時間のスペクトルが何らかの方法で変更される場合、より強い COLA 準拠ウィンドウが推奨されます。

関数 iscola を使用して、弱い COLA 準拠を確認できます。COLA 準拠を確認するために使用される加算数は、ウィンドウの長さとホップ サイズによって決まります。通常、重み付きオーバーラップ加算 (WOLA) には、$${\displaystyle \sum _{m=-\infty }^{\infty }{g}^{a+1}(n-mR)=c}\forall n\in \mathbb{Z},$$ の $$a=1$$、オーバーラップ加算 (OLA) には、$$a=0$$ を使用するのが一般的です。既定では、istft は、オーバーラップ加算方法を実行する前に "合成ウィンドウ" を適用して、WOLA 方法を使用します。

通常、合成ウィンドウは解析ウィンドウと同じです。強い OLA ウィンドウの平方根を取ることで有用な WOLA ウィンドウを構築できます。すべての非負の OLA ウィンドウに対してこの方法を使用できます。たとえば、ルートーハン ウィンドウは WOLA ウィンドウの良い例です。

完全再構成

一般に、入力信号の STFT を計算して逆変換しても、完全再構成は行われません。ISTFT の出力を元の入力信号と可能な限り一致させるには、信号とウィンドウが以下の条件を満たす必要があります。

• 入力サイズ — istft を使用して stft の出力を逆変換して、その結果を入力信号 x と同じ長さにする場合、

  $$k=\frac{N_x-L}{M-L}$$

  の値は整数でなければなりません。方程式の N_x は信号の長さ、M はウィンドウの長さ、L はオーバーラップの長さになります。

• COLA 準拠 — COLA 準拠ウィンドウを使用して、信号の短時間フーリエ変換を変更していないと仮定します。

• パディング — k の値が整数にならないような入力信号の長さの場合、短時間フーリエ変換を計算する前に信号をゼロ パディングします。信号を逆変換した後に、余分なゼロを削除します。

参照

[1] Crochiere, R. E. "A Weighted Overlap-Add Method of Short-Time Fourier Analysis/Synthesis." IEEE Transactions on Acoustics, Speech and Signal Processing. Vol. 28, Number 1, Feb. 1980, pp. 99–102.

[2] Gotzen, A. D., N. Bernardini, and D. Arfib. "Traditional Implementations of a Phase-Vocoder: The Tricks of the Trade." Proceedings of the COST G-6 Conference on Digital Audio Effects (DAFX-00), Verona, Italy, Dec 7–9, 2000.

[3] 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.

[4] Portnoff, M. R. "Time-Frequency Representation of Digital Signals and Systems Based on Short-Time Fourier analysis." IEEE Transactions on Acoustics, Speech and Signal Processing. Vol. 28, Number 1, Feb 1980, pp. 55–69.

[5] Smith, Julius Orion. Spectral Audio Signal Processing. https://ccrma.stanford.edu/~jos/sasp/, online book, 2011 edition, accessed Nov. 2018.

[6] Sharpe, Bruce. Invertibility of Overlap-Add Processing. https://gauss256.github.io/blog/cola.html, accessed July 2019.

拡張機能

tall 配列

メモリの許容量を超えるような多数の行を含む配列を計算します。

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

InputTimeDimension は常に指定し、"downrows" に設定しなければなりません。

詳細については、tall 配列を参照してください。

C/C++ コード生成

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

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

ConjugateSymmetric 引数は、コード生成でサポートされません。

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

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

ConjugateSymmetric 引数は、コード生成でサポートされません。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

この関数は、スレッドベースの環境を完全にサポートします。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

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

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

ConjugateSymmetric が true に設定されていない限り、すべての虚数部がゼロであっても、出力 x は常に複素数です。

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

バージョン履歴

R2019a で導入

参考

オブジェクト

• istftLayer

関数

• dlistft | iscola | pspectrum | spectrogram | stft | stftmag2sig

トピック

• 時間-周波数ギャラリー