フーリエ変換→逆フーリエ変換をした後の波形の再現度が絶望的

Question

###背景
iPhoneのボイスメモで取り込んだピアノの音色をwav変換し、波形表示をしました(navy)。
これをscript.signal.stft/istftして再びmatplotlibで波形を表示させました(aquamarine)。
該当コードと出力は以下の通りです：

import numpy as np
import wave as wave
import scipy.signal as sp
import matplotlib.pyplot as plt
import math

file=wave.open('/Users/***/Desktop/Musica/doremi.wav')

data=file.readframes(file.getnframes())
data=np.frombuffer(data, dtype=np.int16)
duration=file.getnframes()/file.getframerate()
x1_list=[]
for i in range(len(data)):
	x1_list.append(duration*i/len(data))
x1=np.array(x1_list)
y1=np.array(data)

f, t, stft_i=sp.stft(data, fs=file.getframerate(), window='hann', nperseg=512, noverlap=256)
stft_i=10*np.log(np.abs(stft_i))

t, istft_i=sp.istft(stft_i, fs=file.getframerate(), window='hann', nperseg=512, noverlap=256)
x2=t
y2=istft_i

plt.title('Spectrogram of DoReMi')
plt.xlabel('Times(s)')
plt.ylabel('Intensity(Pa)')
plt.plot(x1, y1, color='navy')
plt.plot(x2, y2, color='aquamarine')
plt.show()

###質問①
re-synthesiseした波形の時間が倍になっています。こちらのtというのはx1とは異なり複素数信号とのことですが、何故全体長まで変わってしまうのでしょうか。

###質問②
そもそもintensityがフーリエ変換/逆変換前と全く違います。何が悪かったのでしょう。

素人質問で恐縮ですが、お力添え頂ける箇所が御座いましたら、ご教授お願い申し上げます。

Comments

[toast-uz]

変換と逆変換に間に
stft_i=10*np.log(np.abs(stft_i))
が入っているからではないでしょうか？

Answer 1

変換と逆変換に間に
stft_i=10*np.log(np.abs(stft_i))
が入っているからではないでしょうか？

特にlogにより全体的に音量を潰しているため、逆変換が平板な形になっているように思えます。

追記

上記修正後に時間幅が異なって見えるのは、wavファイルがステレオであるため、dataの前半が1番目のチャネル、後半が2番目のチャネル、となっているためです。dataをもとに変換逆変換をして作ったx2、y2のグラフは、2つのチャネルを時間軸で追記している形になるため、x1、y1のグラフに対して、倍の時間がかかっているように見えてしまいます。

単純に修正するには、x1、y1を計算する前に、以下のようにして1番目のチャネルだけ取り出すとよいでしょう。本来であればチャネルごとに別のグラフにあらわすようにするとよいでしょう。

data = data[:file.getnframes()]

Comments

[jbpb0]

これ、参考になるかと
https://github.com/EnsekiTT/blog/blob/master/20171026_STFT/STFT.ipynb

なお、上記コード中に振幅が amp/10 未満の周波数をカットしているところがありますが、scipy.signal.istftのExamplesで(わざと加えた)信号のノイズを除去する処理と同じことを(なぜか)やっているだけなので、不要です
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.istft.html

[kay_ventris4]

お二方、有難う御座います。

toast-uz様、まず、stft_i=10*np.log(np.abs(stft_i))の箇所を削除したところ、周期が二倍の同じ波形が出現しました。とりあえずx2=t/2で対応したところ、理想の波形を出現させることには成功しましたが、このtの問題は何に起因したものだったのでしょうか。お手数とは存じますが、ご教授賜れれば幸いです。

jbpb0様、「上記コード中に振幅が amp/10 未満の周波数をカットしているところ」と仰いますのは、具体的にどこの部分でしょうか。自分で書いておいて恐縮ですが、ご教授頂けますと嬉しく存じます。

[jbpb0]

「上記コード中に振幅が amp/10 未満の周波数をカットしているところ」とは、私が参考に挙げた
https://github.com/EnsekiTT/blog/blob/master/20171026_STFT/STFT.ipynb
の
Zxx = np.where(np.abs(Zxx) >= amp/10, Zxx, 0)
のことです

kay_ventris4さんのコードのことではありません
分かりにくくてすみませんでした

[kay_ventris4]

jbpb0様、そういうことでしたか。有難う御座います。ご提示くださったものは大変参考になりました。有難う御座います。

[toast-uz]

時間軸がずれている点の修正も追記しました。

[kay_ventris4]

なるほど、np.linspaceを使うといいのですね。
しかしながら、僕が用いたfor文での該当コードの様な記述では誤りということになるのでしょうか？立て続けに何度も申し訳ございません。

[toast-uz]

追記が勘違いしていたので修正しました。計算そのものは私のやり方と質問者様のやり方は変わりありませんでした。

[kay_ventris4]

謎の半分はいつの間にかステレオになっていた片方の音だったんですね。
何度も申し訳ございませんでした。有難うございました。