■追記
ヒートマップを描画するアルゴリズムに大きな改定はありませんが、ヒートマップ描画に使うデータを作成する段階（DFT風演算）で誤りがありました。改訂版は以下に投稿しています。最終的なヒートマップは以下ごご参照ください。
DFT風アルゴリズムバグ対応版でカラスの鳴き声ヒートマップを作る - chakokuのブログ(rev4)
背景：カラスの声をSYNTHで合成したいので周波数を分析する
これまでの取り組み：PCM録音して、Audacityで鳴き声の部分を切り出して、PythonでDFT風に時間ごとの周波数強度を算出(CSV形式で出力)
今日のやること：ヒートマップ（というのか）の形式で、時間ごとの周波数強度をグラフィックで表示
結果：作った
画素が200x200ぐらいで作っていて非常に荒いのだが、、縦軸は周波数、下が1000Hz、上が2000Hz、横軸は時間経過
いろいろ裏付け不十分なまま、とにかく結果を知りたくてなんちゃって作っており、このソース類がどれぐらい正確か？は別の方法で検証必要だが*1、、得られた結果によると、鳴き始めは1300Hzが強くて、途中から1600Hzも含まれる、後半は、1400〰1500Hzが強く発せられるようである。体感的には、鳴き声の最初は高い声から下がるように思えるのだが
周波数取り出しの刻み幅が(Copilotの指摘通り)、100Hz単位で荒いのが問題なのかもしれない。
ソースは以下(WAVは直接読めず、DFT風に周波数強度を演算したCSVを読み込んで、最大値で正規化して赤色の濃さに変換)
#!/usr/bin/python3
import pdb
import csv
import colorsys
from PIL import Image, ImageDraw, ImageFont
CSV_FILE = 'dft_kai.csv'
OUT_PNG_FILE = "h_map.png"
def read_csv(file_name):
rows = []
with open(file_name, newline='', encoding='utf-8') as f:
reader = csv.DictReader(f)
for row in reader:
rows.append(row)
return rows
def get_max_spectrum(dft_data):
max_spectrum = 0.0
header_info = dft_data[0].keys()
for spectrum_in_dic in dft_data[1:]: # for skip header, [1:]
for value_key in header_info:
if value_key == 'nth_sample' or value_key == '':
continue
else:
spectrum = float(spectrum_in_dic[value_key])
if spectrum > max_spectrum:
max_spectrum = spectrum
return max_spectrum
# HSV
HSV_RED = (0.0, 1.0, 1.0)
HSV_PALE_PINK = (0.0, 0.01, 1.0)
HSV_HUE_RED = 0
def hsv_to_rgb255(h, s, v):
r, g, b = colorsys.hsv_to_rgb(h, s, v)
return int(r*255), int(g*255), int(b*255)
RECT_SIZE = 5 # size of heat map rectangle width,height is 5
def draw_heat_map(dft_data, out_file):
header_info = dft_data[0].keys()
#pdb.set_trace()
freq_info = list(header_info)[1:]
# create image 100x100
heat_map_width = 30 + RECT_SIZE * len(dft_data[1:]) # 30 is hanrei for skip header, [1:]
heat_map_height = RECT_SIZE * len(freq_info) #
img = Image.new("RGB", (heat_map_width, heat_map_height), (200, 200, 200))
draw = ImageDraw.Draw(img)
max_spectrum = get_max_spectrum(dft_data)
# draw hanrei
x = 5
y = heat_map_height - RECT_SIZE
font = ImageFont.truetype("/cygdrive/c/Windows/Fonts/meiryo.ttc", 8)
for freq in freq_info:
freq_str = f'{freq}'
draw.text((x, y - 2), freq_str, fill=(0, 0, 0), font=font) # -2 means adhoc adjust ascent / descent
y -= RECT_SIZE
x = 30
y = heat_map_height - RECT_SIZE
for spectrum_in_dic in dft_data[1:]: # for skip header, [1:]
y = heat_map_height - RECT_SIZE
for spectrum in freq_info:
if spectrum == '':
continue
else:
spectrum = float(spectrum_in_dic[spectrum])
saturation = spectrum / max_spectrum
r,g,b = hsv_to_rgb255(HSV_HUE_RED, saturation ,1)
rect = (x, y, x + 5, y + 5)
draw.rectangle(rect, fill=(r, g, b), outline=(200,200,200), width=1)
y -= RECT_SIZE
x += RECT_SIZE
img.save(out_file)
print('output:', out_file)
dft_data = read_csv(CSV_FILE)
draw_heat_map(dft_data, OUT_PNG_FILE)■追記
カラスの声紋？というのか、周波数分析してグラフィック表示しているサイトと比較
CARRION CROW (Corvus corone) Cornelle Noire — wildechoes
（勝手に引用していますが・・・）
ここで掲載されているグラフと比べると最も近い周波数帯は1000Hz～1500Hzなので近いことは近い（上記サイトでは分析している周波数の幅が広く、刻み幅が細かくて高精度なのだが）
■追記
（検証していませんが）だいたい計算できていると想定して、周波数の刻み幅を10Hzにして、調査対象を10Hz～4000Hzまで広げた。この結果は以下(１つの周波数強度を∑で計算する時のサンプリング数は100、サンプリングレート48,000Hzなので、2ms単位で周波数演算していることになる。ウインドウが狭すぎる？？)
（縦軸が周波数、横軸が時間経過、濃い赤色が強く検出された周波数）
窓関数とか使ってないので、縦の周波数方向に滲んでいるのではと思えますが、、これだけ広がっているとSYNTHによる合成が難しい。というか、特徴が見つけられないというか。。 どの周波数帯が強くでているか？は分かるとしても鳴き声の特徴までは見つけられない。Copilotの提案に従って、窓関数(Hann（ハニング）窓)で∑する時のサンプルの両端を正しく処理してから演算してみるか。。
*1:画像の上部に無意味な空間ができている点、凡例の1000Hzが右にずれている点はバグ
