scipy.signal.stftからtorchaudio.transforms.spectrogramへのスペクトログラム計算のコード変換

Question

現在wavファイルを読み込みスペクトログラムに変換するプログラムを書いています。
scipyでプログラムがあるのですが、機械学習を行うためtorchaudioでの実装を試しています。nfftやwindowsizeなど可能な限りパラメータをそろえたのですがスペクトログラムが同じようになりません。返される値のサイズが違うのではと思いD=D/1000を入れてみたところスペクトログラムは同じようにすることができたのですが、なぜ1/1000で同じようになったのかという理屈が分かりません。また、表示の部分でscipyだとstftのときにDの他にもt,ｆが返されますが、torchaudioだとDのみが返されるため縦軸、横軸の表示がうまくいきません。
tに関してはDの時間軸の配列に1/Fsをかければいいというのは分かるのですが、縦軸の周波数成分の配列を用意する方法が分かりません。
わかりにくくなってしまったのですが、
・なぜD/1000にするとうまくいくのか
・scipyでのｆという周波数の配列をtorchaudioもしくはpytorchで作れるのか
以上の2点について教えていただきたいです。
よろしくお願いします。

##use scipy.signal
# read the wav file
ls_org, Fs = sf.read(infilename)
amp = 0.003
nt = min(2**17, 2**(nextpow2(ls_org.size)-1))

ls_org = ls_org[:nt]
ls_org = ls_org / max(ls_org) *0.4
t = np.arange(len(ls_org.T))/Fs
print("Shape of waveform: {}".format(ls_org.shape))
print("Sample rate of waveform: {}".format(Fs))

#compute spectrogram
windowsize = 2048
offset = np.floor(windowsize / 4)

#compute complex spectrogram
f, t, D = signal.stft(ls_org, fs = Fs, nperseg = windowsize , noverlap = windowsize - offset)
(m, n) = D.shape

#Plot signal and spectrogram
fig = pl.figure(1)
pl.plot(ls_org)
pl.xlim([0, len(ls_org)])
pl.title("Input signal", fontsize = 20)
fig = pl.figure(2)
pl.pcolormesh(t, f, np.abs(D), cmap='jet', vmin=0, vmax=0.001)
pl.title('STFT Magnitude')
pl.ylabel('Frequency [Hz]')
pl.xlabel('Time [sec]')

##use torchaudio
#waveread
waveform, Fs = torchaudio.load(infilename)

nt = min(2**17, 2**(nextpow2(len(waveform.T))-1))
waveform = waveform[:, :nt]
waveform = waveform / max(waveform.T) * 0.4
print("Shape of waveform: {}".format(waveform.size()))
print("Sample rate of waveform: {}".format(Fs))
t = np.arange(len(waveform.T))/Fs
#transform spectrogram
nfft = 2048
windowsize = nfft
offset = windowsize // 4
print(offset)
D = torchaudio.transforms.Spectrogram(nfft, windowsize, offset, power=None)(waveform)
D = D / 1000
print("Shape of spectrogram: {}".format(D.size()))
X = torch.view_as_complex(D).clone().detach() #コピー
X = X[0,:,:].numpy()
time = np.arange(257)

#plot
fig = plt.figure(1)
plt.xlabel("time [s]")
plt.ylabel("amplitude")
plt.xlim(0,len(waveform.T)/Fs)
plt.plot(t, waveform.t().numpy())
plt.show()
fig = plt.figure(2)
plt.pcolormesh(time, f, np.abs(X), cmap='jet', vmin=0, vmax=0.001)
plt.xlabel("time [s]")
plt.ylabel("frequency [Hz]")

Answer 1

・なぜD/1000にするとうまくいくのか

1000で割るのではなくて、窓関数の総和で割ります
お書きのコードでは窓関数を指定してないので、scipy.signal.stft()もtorchaudio.transforms.Spectrogram()もデフォルト設定のhann窓が使われてますから、

signal.hann(windowsize).sum()=1023.5

で割ることになります
(ただし、n_fftとwin_lengthが異なる場合にどちらで計算するのかとかは、調べてないので分かりません)

割らないと合わない理由は、scipy.signal.stft()は窓関数の総和で割って正規化する仕様だけど、torchaudio.transforms.Spectrogram()の仕様はそうではないからです
stftに限らず一般のfftも含めて、f特の正規化はこうしなくてはいけないという決まりはないので、ツールによっていろいろです
なので、異なるツールで計算したf特の比較をする際は、それぞれの仕様を調べて適切に換算する必要があります

ちなみに、下記のコードを実行したら分かりますように、numpy.fft.fft()も(デフォルトの設定では)正規化しません

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

xxx = np.arange(0, 2*np.pi*10, 2*np.pi/100)
yyy = np.sin(xxx)

zzz = np.fft.fft(yyy*signal.hann(len(yyy)))
fff = np.fft.fftfreq(len(zzz), d=1/len(zzz))
plt.plot(fff, np.abs(zzz))
plt.xlim(0, 20)
plt.show()
np.max(np.abs(zzz))  # 最大250

fff2, xxx2, zzz2 = signal.stft(yyy, fs=len(yyy))
plt.plot(fff2, np.abs(zzz2[:, int(round(len(xxx2)/2))]))
plt.xlim(0, 20)
plt.show()
np.max(np.abs(zzz2))  # 最大0.44

tに関してはDの時間軸の配列に1/Fsをかければいい

元信号はそうですが、torchaudio.transforms.Spectrogram()の結果は違います
torchaudioを用いたwavファイルのスペクトログラム変換
の回答に書いたように、時間は 0, offset/Fs, 2*offset/Fs... なので、

t2 = np.arange(D.size(1))*offset/Fs

とすれば、scipy.signal.stft()の返り値と同じものを作ることができます

・scipyでのｆという周波数の配列をtorchaudioもしくはpytorchで作れるのか

そちらもやはり
torchaudioを用いたwavファイルのスペクトログラム変換
の回答に書いたように、周波数は 0～Fs/2 なので、

f2 = np.linspace(0, Fs/2, D.size(0))

とすれば、scipy.signal.stft()の返り値と同じものを作ることができます

どうしても torch.* みたいな関数で作るのでないとヤダ、ということでなければ、上記のようにして作れますから、

plt.pcolormesh(t2, f2, np.abs(X), cmap='jet', vmin=0, vmax=0.001)

みたいにできます

Comments

[takumi1114o]

回答ありがとうございます。
正規化するかどうかが違ったのですね！
本当にわかりやすい説明ありがとうございました！
jbpb0さんのおかげで理解しながらうまく動かせました