PythonとRを用いたFFTについて

Question

前提・実現したいこと

wavファイルにFFTをかけてスペクトルを表示するコードを
PythonとRの２パターンで書きました。

発生している問題・エラーメッセージ

Pythonでは以下のスペクトルが表示されました。

Rでは以下のスペクトルが表示されました。

どちらも同じような形のスペクトルが表示されているのですが、
ピークが出ている周波数に倍くらいズレが出ています。

Pythonでは140Hzあたりに一番大きなピークがあるのですが、
Rでは280Hzあたりに一番大きなピークが出ています。

どちらのコードが間違っているかわからず、困っている状況です。

該当のソースコード

import wave
import numpy as np
import scipy.fftpack
import matplotlib.pyplot as plt

wf = wave.open("~~~.wav" , "r" )
fs = wf.getframerate()  # サンプリング周波数
x = wf.readframes(wf.getnframes())
x = np.frombuffer(x, dtype= "int16") / 32768.0  # -1 - +1に正規化
wf.close()

start = 0  # サンプリングする開始位置
N = 1550000 # FFTのサンプル数

X = np.fft.fft(x[start:start+N])  # FFT

freqList = np.fft.fftfreq(N, d=1.0/fs)  # 周波数軸の値を計算

plt.title("FFT-wheeze")
plt.plot(freqList,abs(X))
plt.axis([0,600,0,max(abs(X))+500])

plt.show()

install.packages("tuneR")
library(tuneR)

# ファイルの読み込み
wavobj <- readWave("~~~.wav")
x <- wavobj@left         # 左チャンネルだけを使う（wavobj@.Dataを使う場合も）
fs <- [e-mail]   # 標本化周波数
nbits <- wavobj@bit      # 量子化ビット数

x <- x[1:(fs*5)]   # 先頭の5秒間だけを分析対象として抜き出す

y <- fft(x)

y.tmp <- Mod(y)   # 複素数の絶対値を計算

# 振幅スペクトル（amplitude spectrum）の計算
y.tmp <- Mod(y)   # 複素数の絶対値を計算
y.ampspec <- y.tmp[1:(length(y)/2+1)]   # DC〜ナイキストの範囲
y.ampspec[2:(length(y)/2)] <- y.ampspec[2:(length(y)/2)] * 2

# 対応する周波数の計算（DC〜ナイキストの範囲をビン数で均等割）
f <- seq(from=0, to=fs/2, length=length(y)/2+1)

# グラフ作成
plot(f, y.ampspec, type="h", xlab="Frequency (Hz)", ylab="Amplitude Spectrum", xlim=c(0, 600))
grid()

Comments

[tachikoma]

wavはLRの2チャンネルですか？

Answer 1

Pythonのwavモジュールをそのまま使うなら、LRの分割を自分でやる必要があります。おそらく左右にほとんど同じ信号が入っているため、分割前だと周波数が半分になったんだと思います。

wf = wave.open("~~~.wav" , "r" )
fs = wf.getframerate()  # サンプリング周波数
nchannels = wf.getnchannels()
x = wf.readframes(wf.getnframes()*nchannels)
wf.close()

# bytes arrayからnumpy.arrayに変換
x = np.frombuffer(x, dtype= "int16") / 32768.0

# stereo音源の場合、左右に分割
if nchannles == 2:
    Lchannle = x[::2]
    Rchannle = x[1::2]

Comments

[k.n.]

解決しました。ありがとうございます。
1つわからないのですが、
x = wf.readframes(wf.getnframes()*nchannels)でnchannelsをかけているのはなぜですか。

[tachikoma]

nchannelが必要なのは、readframesの実装がbit depthは考慮するけど、チャネル数を考慮してくれないためです。monoだろうがstereoだろうが同じだけだけのデータ数を読み込みます。(wavのフォーマットにもよりますが)stereoだとLRLR...の順にデータが並んでおり、monoだとLLLL....と並んいるので、4frames読み込む場合、stereoだとLRLRを読み込むのに対して、monoだとLLLLを読み込みます。このため、monoと同じだけの時間のデータをstereoから読み込むには読み込むframes x nchannelsが必要なわけです。

[tachikoma]

補足ですが、音声ファイルの読み込みはsoundfileを使うと便利ですよ。ご参考までに。

[k.n.]

丁寧にありがとうございます。理解できました。
soundfileコード簡潔でいいですね。参考にします。