FFTとエネルギーについての考察 - TadaoYamaokaの開発日記

前回の日記で、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