PythonとRを用いたFFTについて
解決済
回答 1
投稿
- 評価
- クリップ 0
- VIEW 2,990
前提・実現したいこと
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 <- wavobj@samp.rate # 標本化周波数
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()-
気になる質問をクリップする
クリップした質問は、後からいつでもマイページで確認できます。
またクリップした質問に回答があった際、通知やメールを受け取ることができます。
クリップを取り消します
-
良い質問の評価を上げる
以下のような質問は評価を上げましょう
- 質問内容が明確
- 自分も答えを知りたい
- 質問者以外のユーザにも役立つ
評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。
質問の評価を上げたことを取り消します
-
評価を下げられる数の上限に達しました
評価を下げることができません
- 1日5回まで評価を下げられます
- 1日に1ユーザに対して2回まで評価を下げられます
質問の評価を下げる
teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。
- プログラミングに関係のない質問
- やってほしいことだけを記載した丸投げの質問
- 問題・課題が含まれていない質問
- 意図的に内容が抹消された質問
- 過去に投稿した質問と同じ内容の質問
- 広告と受け取られるような投稿
評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。
質問の評価を下げたことを取り消します
この機能は開放されていません
評価を下げる条件を満たしてません
質問の評価を下げる機能の利用条件
この機能を利用するためには、以下の事項を行う必要があります。
- 質問回答など一定の行動
-
メールアドレスの認証
メールアドレスの認証
-
質問評価に関するヘルプページの閲覧
質問評価に関するヘルプページの閲覧
checkベストアンサー
+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]投稿
-
回答の評価を上げる
以下のような回答は評価を上げましょう
- 正しい回答
- わかりやすい回答
- ためになる回答
評価が高い回答ほどページの上位に表示されます。
-
回答の評価を下げる
下記のような回答は推奨されていません。
- 間違っている回答
- 質問の回答になっていない投稿
- スパムや攻撃的な表現を用いた投稿
評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。
15分調べてもわからないことは、teratailで質問しよう！
- ただいまの回答率 88.11%
- 質問をまとめることで、思考を整理して素早く解決
- テンプレート機能で、簡単に質問をまとめられる
質問への追記・修正、ベストアンサー選択の依頼
tachikoma
2018/07/07 10:07
wavはLRの2チャンネルですか？