読者です 読者をやめる 読者になる 読者になる

TadaoYamaokaの日記

山岡忠夫 Home で公開しているプログラムの開発ネタを中心に書いていきます。

FFTとエネルギーについての考察

前回の日記で、FFTのピークの値は周波数分解能の整数倍とのずれによって値が変わるため単純には比較できないことを説明した。

ここでは、周波数全体のエネルギーについては、周波数分解能の整数倍とのずれによらず一定となることについて考察する。

エネルギーが一定になることについてはChachayさん日記で言及されていて気づきました。
chachay.hatenablog.com


解析元の音声データのエネルギーは、振幅の2乗の総和となるが、パーセバルの定理によりFFTの値の2乗(パワースペクトル)の総和で計算できる。

\displaystyle \sum_{n=0}^{N-1} |y_n|^2 = \frac{1}{N} \sum_{n=0}^{N-1} |Y_n|^2


440Hzの正弦波と435Hzの正弦波のエネルギーを比較すると以下の通りとなる。

import numpy as np
import matplotlib.pyplot as plt
fs = 44100
size = 4096
t = np.arange(0, size) / fs
frq = np.fft.fftfreq(size, 1/fs)

#440Hz
f = 440
y = np.sin(2 * np.pi * f * t)
Y = np.fft.fft(y)
P = abs(Y)**2
sum(P)/size

#435Hz
f2 = 435
y2 = np.sin(2 * np.pi * f2 * t)
Y2 = np.fft.fft(y2)
P2 = abs(Y2)**2
sum(P2)/size
440Hz 2051.68
435Hz 2051.62

ほぼ同じ値となっている。
裾の広がり方は変わるが、パワースペクトルを合計すると一致することが確かめられた。


音声データを解析する場合、基本周波数と倍音が混じっているため、基本周波数と倍音のエネルギーを比較したい場合に全周波数の和を計算するわけにはいかない。
しかし、周辺の周波数のパワースペクトルの和で、ある程度比較可能な値になるのではないかと思う。
以下の通り検証してみた。


まず、440Hzと435Hzの正弦波のパワースペクトルのピークを比較すると以下の通りとなる。

#440Hz
np.max(P)
np.max(np.log(P))

#435Hz
np.max(P2)
np.max(np.log(P2))
周波数 パワースペクトル(ピーク) logスケール
440Hz 3964274.05 15.19
435Hz 2401832.86 14.69

1.65倍の差がある。


周辺3つのパワースペクトルの和のピークを比較すると以下の通りとなる。

#440Hz
P_around = np.convolve(P, np.ones(3), 'valid')
np.max(P_around)
np.max(np.log(P_around))

#435Hz
P2_around = np.convolve(P2, np.ones(3), 'valid')
np.max(P2_around)
np.max(np.log(P2_around))

※np.convolveは畳み込み。値が1のベクトルとの畳み込みで和が計算できる。

周波数 周辺のパワースペクトルの和(ピーク) logスケール
440Hz 4111315.81 15.22
435Hz 3674791.85 15.11


差は1.11倍となり、パワースペクトルのピークに比べて小さくなっている。

グラフにすると以下のようになる。

f:id:TadaoYamaoka:20160829223017p:plain
※logスケール


昨日の日記で解析した435Hzと基本周波数と2倍音(振幅95%)、3倍音(振幅90%)を含むデータを周辺のパワースペクトルの和でグラフにすると以下のようになる。

f = 435
y = np.sin(2 * np.pi * f * t) + np.sin(2 * np.pi * 2 * f * t) * 0.95 + np.sin(2 * np.pi * 3 * f * t) * 0.9
Y = np.fft.fft(y)
P = abs(Y)**2
P_around = np.convolve(P, np.ones(3), 'valid')
plt.plot(frq[1:size-1], P_around)
plt.axis([0, 2000, 0, 4000000])
plt.show()

f:id:TadaoYamaoka:20160829223918p:plain

FFTの値では、2倍音のピーク一番高くなっていたが、基本周波数が一番高くなっている。続いて2倍音、3倍音と元のデータの振幅と同じ順になっている。


FFTのピークを比較する場合は、周辺のパワースペクトルの和(もしくは移動平均)を比較することで、ある程度周波数分解能の整数倍とのずれの影響を抑えることができそうである。
実際の音声データで試した場合にうまく機能するかは今後試してみたい。