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

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

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

R言語は、「S言語」をオープンソースとして実装なおした、統計解析向けのプログラミング言語です。 計算がとても速くグラフィックも充実しているため、数値計算に向いています。 文法的には、統計解析部分はS言語を参考にしており、データ処理部分はSchemeの影響を受けています。 世界中の専門家が開発に関わり、日々新しい手法やアルゴリズムが追加されています。

Python

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

Q&A

解決済

1回答

5773閲覧

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

k.n.

総合スコア9

R

R言語は、「S言語」をオープンソースとして実装なおした、統計解析向けのプログラミング言語です。 計算がとても速くグラフィックも充実しているため、数値計算に向いています。 文法的には、統計解析部分はS言語を参考にしており、データ処理部分はSchemeの影響を受けています。 世界中の専門家が開発に関わり、日々新しい手法やアルゴリズムが追加されています。

Python

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

0グッド

0クリップ

投稿2018/07/06 18:17

前提・実現したいこと

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

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

Pythonでは以下のスペクトルが表示されました。
イメージ説明
Rでは以下のスペクトルが表示されました。
イメージ説明

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

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

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

該当のソースコード

Python

1import wave 2import numpy as np 3import scipy.fftpack 4import matplotlib.pyplot as plt 5 6wf = wave.open("~~~.wav" , "r" ) 7fs = wf.getframerate() # サンプリング周波数 8x = wf.readframes(wf.getnframes()) 9x = np.frombuffer(x, dtype= "int16") / 32768.0 # -1 - +1に正規化 10wf.close() 11 12start = 0 # サンプリングする開始位置 13N = 1550000 # FFTのサンプル数 14 15X = np.fft.fft(x[start:start+N]) # FFT 16 17freqList = np.fft.fftfreq(N, d=1.0/fs) # 周波数軸の値を計算 18 19plt.title("FFT-wheeze") 20plt.plot(freqList,abs(X)) 21plt.axis([0,600,0,max(abs(X))+500]) 22 23plt.show()

R

1install.packages("tuneR") 2library(tuneR) 3 4# ファイルの読み込み 5wavobj <- readWave("~~~.wav") 6x <- wavobj@left # 左チャンネルだけを使う(wavobj@.Dataを使う場合も) 7fs <- wavobj@samp.rate # 標本化周波数 8nbits <- wavobj@bit # 量子化ビット数 9 10x <- x[1:(fs*5)] # 先頭の5秒間だけを分析対象として抜き出す 11 12y <- fft(x) 13 14y.tmp <- Mod(y) # 複素数の絶対値を計算 15 16# 振幅スペクトル(amplitude spectrum)の計算 17y.tmp <- Mod(y) # 複素数の絶対値を計算 18y.ampspec <- y.tmp[1:(length(y)/2+1)] # DC〜ナイキストの範囲 19y.ampspec[2:(length(y)/2)] <- y.ampspec[2:(length(y)/2)] * 2 20 21# 対応する周波数の計算(DC〜ナイキストの範囲をビン数で均等割) 22f <- seq(from=0, to=fs/2, length=length(y)/2+1) 23 24# グラフ作成 25plot(f, y.ampspec, type="h", xlab="Frequency (Hz)", ylab="Amplitude Spectrum", xlim=c(0, 600)) 26grid()

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

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

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

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

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

tachikoma

2018/07/07 01:07

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

回答1

0

ベストアンサー

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

python

1wf = wave.open("~~~.wav" , "r" ) 2fs = wf.getframerate() # サンプリング周波数 3nchannels = wf.getnchannels() 4x = wf.readframes(wf.getnframes()*nchannels) 5wf.close() 6 7# bytes arrayからnumpy.arrayに変換 8x = np.frombuffer(x, dtype= "int16") / 32768.0 9 10# stereo音源の場合、左右に分割 11if nchannles == 2: 12 Lchannle = x[::2] 13 Rchannle = x[1::2] 14

投稿2018/07/07 01:19

編集2018/07/07 11:06
tachikoma

総合スコア3601

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

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

k.n.

2018/07/07 06:39

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

2018/07/07 11:16

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

2018/07/07 11:17

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

2018/07/09 07:01

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問