質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
PyTorch

PyTorchは、オープンソースのPython向けの機械学習ライブラリ。Facebookの人工知能研究グループが開発を主導しています。強力なGPUサポートを備えたテンソル計算、テープベースの自動微分による柔軟なニューラルネットワークの記述が可能です。

Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

コードレビュー

コードレビューは、ソフトウェア開発の一工程で、 ソースコードの検査を行い、開発工程で見過ごされた誤りを検出する事で、 ソフトウェア品質を高めるためのものです。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

1回答

1614閲覧

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

takumi1114o

総合スコア2

PyTorch

PyTorchは、オープンソースのPython向けの機械学習ライブラリ。Facebookの人工知能研究グループが開発を主導しています。強力なGPUサポートを備えたテンソル計算、テープベースの自動微分による柔軟なニューラルネットワークの記述が可能です。

Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

コードレビュー

コードレビューは、ソフトウェア開発の一工程で、 ソースコードの検査を行い、開発工程で見過ごされた誤りを検出する事で、 ソフトウェア品質を高めるためのものです。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2021/01/02 15:29

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

Python3

1##use scipy.signal 2# read the wav file 3ls_org, Fs = sf.read(infilename) 4amp = 0.003 5nt = min(2**17, 2**(nextpow2(ls_org.size)-1)) 6 7ls_org = ls_org[:nt] 8ls_org = ls_org / max(ls_org) *0.4 9t = np.arange(len(ls_org.T))/Fs 10print("Shape of waveform: {}".format(ls_org.shape)) 11print("Sample rate of waveform: {}".format(Fs)) 12 13#compute spectrogram 14windowsize = 2048 15offset = np.floor(windowsize / 4) 16 17#compute complex spectrogram 18f, t, D = signal.stft(ls_org, fs = Fs, nperseg = windowsize , noverlap = windowsize - offset) 19(m, n) = D.shape 20 21#Plot signal and spectrogram 22fig = pl.figure(1) 23pl.plot(ls_org) 24pl.xlim([0, len(ls_org)]) 25pl.title("Input signal", fontsize = 20) 26fig = pl.figure(2) 27pl.pcolormesh(t, f, np.abs(D), cmap='jet', vmin=0, vmax=0.001) 28pl.title('STFT Magnitude') 29pl.ylabel('Frequency [Hz]') 30pl.xlabel('Time [sec]') 31 32##use torchaudio 33#waveread 34waveform, Fs = torchaudio.load(infilename) 35 36nt = min(2**17, 2**(nextpow2(len(waveform.T))-1)) 37waveform = waveform[:, :nt] 38waveform = waveform / max(waveform.T) * 0.4 39print("Shape of waveform: {}".format(waveform.size())) 40print("Sample rate of waveform: {}".format(Fs)) 41t = np.arange(len(waveform.T))/Fs 42#transform spectrogram 43nfft = 2048 44windowsize = nfft 45offset = windowsize // 4 46print(offset) 47D = torchaudio.transforms.Spectrogram(nfft, windowsize, offset, power=None)(waveform) 48D = D / 1000 49print("Shape of spectrogram: {}".format(D.size())) 50X = torch.view_as_complex(D).clone().detach() #コピー 51X = X[0,:,:].numpy() 52time = np.arange(257) 53 54#plot 55fig = plt.figure(1) 56plt.xlabel("time [s]") 57plt.ylabel("amplitude") 58plt.xlim(0,len(waveform.T)/Fs) 59plt.plot(t, waveform.t().numpy()) 60plt.show() 61fig = plt.figure(2) 62plt.pcolormesh(time, f, np.abs(X), cmap='jet', vmin=0, vmax=0.001) 63plt.xlabel("time [s]") 64plt.ylabel("frequency [Hz]")

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

0

ベストアンサー

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

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

python

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

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

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

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

python

1import numpy as np 2from scipy import signal 3import matplotlib.pyplot as plt 4 5xxx = np.arange(0, 2*np.pi*10, 2*np.pi/100) 6yyy = np.sin(xxx) 7 8zzz = np.fft.fft(yyy*signal.hann(len(yyy))) 9fff = np.fft.fftfreq(len(zzz), d=1/len(zzz)) 10plt.plot(fff, np.abs(zzz)) 11plt.xlim(0, 20) 12plt.show() 13np.max(np.abs(zzz)) # 最大250 14 15fff2, xxx2, zzz2 = signal.stft(yyy, fs=len(yyy)) 16plt.plot(fff2, np.abs(zzz2[:, int(round(len(xxx2)/2))])) 17plt.xlim(0, 20) 18plt.show() 19np.max(np.abs(zzz2)) # 最大0.44

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

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

python

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

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

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

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

python

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

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

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

python

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

みたいにできます

投稿2021/01/04 07:00

jbpb0

総合スコア7651

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

takumi1114o

2021/01/05 05:47

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問