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

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

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

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

Q&A

解決済

1回答

742閲覧

Pythonで音声解析

HiruLow

総合スコア55

Python

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

0グッド

0クリップ

投稿2018/03/11 22:55

Pythonを使ってWavデータからピッチを取得する関数を描いているのですが、以下の関数の最後に出る p とゆう変数にWavデータを解析した結果が出る予定だったのですが、該当する値が出てきません。

試しに100Hzを鳴らした音声を検索した際の結果は10163.6634004となってしまいました。
アドバイスいただけたらと思います。

#coding:utf-8 import wave import numpy as np import matplotlib as mpl mpl.use('Agg') import matplotlib.pyplot as plt wf = wave.open("./100hz.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 = 512 # FFTのサンプル数 SHIFT = 128 # 窓関数をずらすサンプル数 hammingWindow = np.hamming(N) freqList = np.fft.fftfreq(N, d=1.0/fs) # 周波数軸の値を計算:Hz x値 def get_hz(x): maxV = 0 maxN = 0 i = 0 for v in x: if v > maxV: if v < 1: maxV = v maxN = i i+=1 return maxN def search(): global start windowedData = hammingWindow * x[start:start+N] # 切り出した波形データ(窓関数あり) X = np.fft.fft(windowedData) # FFT amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X] # 振幅スペクトル:Y値 time=start+N#Time:t値 hznum = 0 hznum = get_hz(amplitudeSpectrum) tmp = freqList[hznum] if tmp < 0: tmp += tmp tmp *= -1 freq = tmp if hznum > 0: if hznum < len(amplitudeSpectrum)-1: tmp1 = freqList[hznum - 1] tmp2 = freqList[hznum + 1] tmp3 = freqList[hznum] if tmp1 < 0: tmp1 += tmp1 tmp1 *= -1 if tmp2 < 0: tmp2 += tmp2 tmp2 *= -1 if tmp3 < 0: tmp3 += tmp3 tmp3 *= -1 dL = tmp1 / tmp3 dR = tmp2 / tmp3 freq += 0.5 * (dR * dR - dL * dL) p = freq * (len(amplitudeSpectrum)/2) / len(amplitudeSpectrum) start += SHIFT # 窓関数をかける範囲をずらす if start + N > len(x): () else: () search()

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

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

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

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

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

退会済みユーザー

退会済みユーザー

2018/03/11 23:37

振幅最大の周波数を知りたい、ということでしょうか?
HiruLow

2018/03/12 01:22

そうゆう事になります。音声はアカペラの歌声等で、その中から発生している音階を取得する為のHzを算出させたいです
退会済みユーザー

退会済みユーザー

2018/03/12 01:27

コードをまだ読んでないですが、librosaを使ったら楽をできるかもしれません。 https://librosa.github.io/librosa/index.html 取り急ぎ参考情報までに。
HiruLow

2018/03/12 02:34

librosaについても検討したのですが、描画処理はすぐれていますが情報を取得するにはあまりむいていませんでした。
guest

回答1

0

ベストアンサー

基本アイデアはdkato0077さんの質問コメントの通りだろうと思いますが、最大スペクトルの両隣りも見ているところから基本周波数(fs/N)の倍数でない周波数の推定も狙っておられるようですね。

###最大のスペクトルを調べる範囲
FFTの扱いをわざわざ面倒にしているような気がします。FFTの結果のスペクトルのうちナイキスト周波数より大きな部分の成分には数学的な意味はあるものの取り扱いが厄介なだけな気がするのでとりあえず無視してはどうでしょう?ナイキスト周波数(=サンプリング周波数の半分)以上の周波数の位相はサンプリング区間ごとにπ以上変化することになります。DFTの考え方を大雑把に捉えればサンプルごとに位相変化がπ未満(理想的にはそれよりずっと小さい)を期待する計算方法と言えると思います。

具体的にはfreqListの前半分のところだけ処理対象にした方がよいと思います。そうすればナイキスト周波数未満の周波数のみになりfreqListの要素も全て正の周波数になってます。

###基本周波数からのずれの推定
元の周波数を推定する計算をしておられますが、計算式が間違っているような気がします。

0.5 * (dR * dR - dL * dL)という式を自分が見たことあるのは以下のページですけど

https://answers.unity.com/questions/157940/getoutputdata-and-getspectrumdata-they-represent-t.html

この計算の意味は「スペクトル最大の周波数の補正」ではなく「スペクトル最大の周波数のDFT結果のインデックスの補正」だと思います。ちなみに自分もこの数式の厳密な説明はできず、3点を通過する曲線を想定してその最大値の場所を推定する式なんかなぁと漠然と推測する程度なので、本当に目的に合ってるかはよく分かりません...


質問者さんのコードに基づきトライしてみました。
wavファイルを使っておられますが、論理を確認する段階ではもっと単純にデータを自動生成するなど工夫しておいた方がデバッグがしやすいと思います。以下ではサンプリング周波数44.1KHz、解析対象の周波数を100Hzとしてやってみてます。

このようにしておけば、再現が容易ですので、閲覧する立場としても問題点の核心に迫り易い気がします。

Python

1import numpy as np 2import matplotlib.pyplot as plt 3 4fs = 44100 5total_sec = 1 6target_freq = 100 7t = np.arange(0, total_sec, 1 / fs) 8x = np.sin((2 * np.pi * target_freq) * t) 9 10start = 0 # サンプリングする開始位置 11N = 512 # FFTのサンプル数 12 13hammingWindow = np.hamming(N) 14freqList = np.fft.fftfreq(N, d=1.0 / fs)[:N // 2] # 周波数軸の値を計算:Hz x値 15 16 17def search(): 18 windowedData = hammingWindow * x[start:start+N] # 切り出した波形データ(窓関数あり) 19 X = np.fft.fft(windowedData)[:N // 2] # FFT : 結果の後半は捨てる 20 amplitudeSpectrum = np.abs(X) # 振幅スペクトル:Y値 21 hznum = np.argmax(amplitudeSpectrum) 22 23 if 0 < hznum < len(freqList) - 1: 24 spectramL = amplitudeSpectrum[hznum - 1] 25 spectramM = amplitudeSpectrum[hznum] 26 spectramR = amplitudeSpectrum[hznum + 1] 27 dL = spectramL / spectramM 28 dR = spectramR / spectramM 29 hznum += 0.5 * (dR * dR - dL * dL) 30 31 p = hznum * fs / N 32 print('p=', p) 33 return p, freqList, X, x[start:start+N], windowedData 34 35 36p, xx, yy, rss, wss = search() 37 38g1 = plt.subplot("311") 39g1.set_title('wave (raw)') 40g1.plot(t[:len(rss)], rss) 41 42g2 = plt.subplot("312") 43g2.set_title('wave (by hamming window)') 44g2.plot(t[:len(rss)], wss) 45g3 = plt.subplot("313") 46 47g3.set_title('spectrum') 48si_min, si_max = 0, 5 49g3.plot(xx[si_min:si_max], np.abs(yy[si_min:si_max]), marker='o') # 注目する部分だけ拡大 50g3.scatter([p], [1], color='red', marker='*') # 周波数推定値をプロット 51 52plt.tight_layout() 53plt.show()

イメージ説明

赤い'*'で示した点が推定周波数です。サンプリング周波数に対して対象周波数が低すぎて誤差が大きくなっているせいか、はたまたなにかバグがあるのか分かりませんが100Hzピッタリにはなってないですね。

投稿2018/03/12 06:54

編集2018/03/12 23:36
KSwordOfHaste

総合スコア18392

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

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

退会済みユーザー

退会済みユーザー

2018/03/12 07:07

44.1kのサンプリング周波数での100Hzは一周忌が441点なので、N=512じゃギリギリ一周期分しか入らないので欲しい点の周りがまばらになり過ぎてしまいますぉ。N=4096とかにまで広げれば少しはマシになると思います。
KSwordOfHaste

2018/03/12 07:10

そうですね。wavのサンプリング周波数がいくつかがわからなかったので敢て「誤差はこんなものかも知れない。しかしそれでも生のFFT結果よりはずっとナイスな推定ができる」といった感じのコメントにしてみました。
退会済みユーザー

退会済みユーザー

2018/03/12 07:13

なるほどなるほど。
KSwordOfHaste

2018/03/12 08:09

質問者さんへ: サンプリング周波数はそのままにFFT区間を変えてみると精度があがる様子が観察できます。♩=120のテンポでFFT区間が四分音符の何倍かを添えてみますと 512 : 96.20 (0.02♩) 2048: 101.32 (0.09♩) 4096 : 99.26 (0.18♩) 8192 : 100.25 (0.37♩=>これくらいになると8文音符とか16分音符の解析が厳しい・・・) 逆にFFT区間を固定してサンプリング周波数を変えてみると 22050 : 96.94 11025 : 101.32 4000 : 100.45 人間の歌声を解析するならサンプリング周波数はあまり高くしない方がよさそうです。
HiruLow

2018/03/12 08:14

とても有力な情報ありがとうございます。 各サイトを見ても日本語でなかったり、式だけが漠然と乗っていたりで、なかなか勉強しようにも理解できないでいました。とてもわかりやすく解析の参考になります。 ありがとうございます。
KSwordOfHaste

2018/03/12 08:17

バリトン(100Hz?)~ソプラノ(1500?)ぐらいがカバーできるのを狙ってると思います。FFT区間とサンプリング周波数を調整しつつ誤差を調べてみるとよいのではないでしょうか。音声を一々用意するのは面倒だと思うので最初は自動生成した波形で当たりを付けるのがよいかも知れません。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問