■追記
本書き込みではDFT風に周波数の強さを計算していますが、アルゴリズムが間違っていることが分かっています。近いうちに訂正します(専門書による理解ではなく、Copilotと相談して得られた情報です)
- 間違っている点１：切り取りサイズが1000サンプルの時、計算できる周波数の最低と最高が決まるはずですがそれを考慮せずに計算できないはずの周波数まで計算している（サンプリングレート48KHzから取り出した1000サンプル分なので1000サンプルにきっちり入る1周期の周波数は48Hz、調べられる最低周波数は48Hz、最高周波数は24KHz(ナイキスト周波数)）
- 間違っている点２：切り取りサイズが1000サンプルの時、計算できる周波数は、1000サンプルを1周期とした周波数(48Hz)に対する整数倍しか調べられないはずなのに、48Hzの整数倍を無視していろんな周波数を計算している
- 機能不足１：調査対象の周波数が上記48Hzの整数倍でない場合、sin/cosの波がデータの両端ブチ切りになってスペクトルリークが発生する。48Hzの整数倍以外で周波数を調査するには、HANN窓を掛け算して両端が滑らかに０になるように補正する
- 機能不足２：DFTによる信号の強さはサンプル数に比例する。正規化するにはサンプル数(N=1000)で割り算する
背景：デジタル録音されたカラスの鳴き声の周波数成分を時間ごとに取り出したい
アプローチ：本来はPythonのDFTのライブラリで計算すべきだが、、DFT風に自力で計算してみる
結果：いろいろ試して（恣意的に良いパラメータで計算すると）、発話時間によって、1300Hz～1600Hzが強くでているようである
詳細：
フーリエ級数の考え方をベースに、調べたい周波数ごとにsin/cosを作って鳴き声データに離散的に掛け算して∑する。DFT風のプログラムは以下。高速化とか考えていなくて、力業で計算。
#!/usr/bin/python3
#
# dft like
#
import math
IN_FILE = 'sound.csv'
SAMPLING_FREQ = 48000
def get_sin_cos(freq, sample, nth):
val_sin = math.sin(2 * math.pi * freq / sample * nth)
val_cos = math.cos(2 * math.pi * freq / sample * nth)
return(val_sin, val_cos)
nth_val = 0
freq_inc= {}
for freq in range(1000, 2001, 100):
freq_inc[freq]=0
print('nth_sample, ', end='')
for freq in range(1000, 2001, 100):
print(freq, end=',')
print('')
with open(IN_FILE,'r') as f:
for line in f:
l_ch, r_ch = line.split(',')
l_ch = int(l_ch)
r_ch = int(r_ch)
for freq in range(1000, 2001, 100):
sin, cos = get_sin_cos(freq, SAMPLING_FREQ, nth_val)
freq_inc[freq] += l_ch * sin + l_ch * cos
if nth_val % 1000 == 0:
print(nth_val,end=', ')
for freq in range(1000, 2001, 100):
print(freq_inc[freq], end=', ')
freq_inc[freq]=0
print('')
nth_val += 1 上記はCSV形式データなので、Excelに読み込ませてグラフ化する。上記以外の条件でも計算できるが、窓関数とか使っていないせいかあまりはっきりとした傾向が出せない。見て分かりやすい条件で得られた結果は以下
鳴き始めてから系列４（1300Hz）が強くなって、中盤は系列５（1400Hz）と系列６（1500Hz）が強くなると判断できる。スマフォのFFTで計測した時は、1400Hzが強く波形出ていたので、DFT風解析でも同じ結果になっている。
今回はCSVに出力したものをDFT風に周波数ごとの強さを調べたのだが、これをPythonだけでヒートマップのような点の大きさとか色の濃さで表現したい。
■追記
上記コードを生成AIに合っているか？と聞いてみた。DFTを使う考え方はあっているが、cosは実数部で、sinは虚数部であり、これらの大きさの計算は、直接加算してはだめで、2乗して加算してSQRT取れと言ってきた*1。やはり理解不十分なのであった。そのようにした結果は以下
改訂版ソースは以下。生成AIからは、音声入力Lだけ使うのではなくて、L/Rの平均を取れとか、窓関数使えとか、周波数の刻み幅が大きすぎるとか、いろいろ指摘を受けた。さすが。。
#!/usr/bin/python3
#
# dft like v0.1
#
# modify: sqrt(real_number ** 2 + im_number ** 2)
import math
IN_FILE = 'sound.csv'
SAMPLING_FREQ = 48000
def get_sin_cos(freq, sample, nth):
val_sin = math.sin(2 * math.pi * freq / sample * nth)
val_cos = math.cos(2 * math.pi * freq / sample * nth)
return(val_sin, val_cos)
nth_val = 0
freq_rel = {}
freq_im = {}
for freq in range(1000, 2001, 100):
freq_rel[freq]=0
freq_im[freq]=0
print('nth_sample, ', end='')
for freq in range(1000, 2001, 100):
print(freq, end=',')
print('')
with open(IN_FILE,'r') as f:
for line in f:
l_ch, r_ch = line.split(',')
l_ch = int(l_ch)
r_ch = int(r_ch)
for freq in range(1000, 2001, 100):
sin, cos = get_sin_cos(freq, SAMPLING_FREQ, nth_val)
freq_rel[freq] += l_ch * cos
freq_im[freq] += l_ch * sin
if nth_val % 1000 == 0:
print(nth_val,end=', ')
for freq in range(1000, 2001, 100):
amp = math.sqrt(freq_rel[freq]**2 + freq_im[freq]**2)
print(amp, end=', ')
freq_rel[freq]=0
freq_im[freq]=0
print('')
nth_val += 1*1:フーリエ変換の考え方は、直交する新しい座標軸を導入し値を写像することであり、元の値にsin/cosで掛け算して∑する(ここは離散的な計算で代用)ことで、新しい座標軸の各軸における値が得られて、プロットすることができる。周波数の類似度（大きさ）とは、プロットされた位置の原点からの距離なので、２つの座標の2乗のSQRTであると。。
