カラスの声を周波数分析ー＞ヒートマップ作成（改良版）

背景：カラスの声をSYNTHで合成したい。そのためまずは声の特徴を周波数の観点で見える化
課題：DFT風演算する際、オーバーラップせずにN=1000で取り出していた
取り組み：Copilotの提案に従い、50%ずつオーバーラップさせながら、N=1000のWindowサイズでDFT風変換して周波数の強さを算出
結果：横方向に2倍の細かさでヒートマップを作成した。時間経過における周波数の変化が前回より細かく表現可能になった。

図中、上がカラスの声(LchとRch)、下が周波数ごとの強さを表すヒートマップ

ソースは以下(周波数強度を出すところまで。ヒートマップの作成は前回のソースのまま流用)

#!/usr/bin/python3

#
# DFT like Spectrum v0.01 (2026/3/4)
#

import pdb
import struct
import math


IN_FILE='260227_0038.wav'
OUT_FILE='dft_out.csv'

WINDOW_SIZE = 1_000   #   lowest frequency is 48Hz

SAMPLING_FREQUENCY = 48_000

FREQUENCY_RESOLUTION = int(SAMPLING_FREQUENCY / WINDOW_SIZE)

MIN_ANALYZE_FREQ = FREQUENCY_RESOLUTION
MAX_ANALYZE_FREQ = int(2000 / FREQUENCY_RESOLUTION) * FREQUENCY_RESOLUTION  # 2000 is depended on user's spec


def read_WAV_headers(f):

    # read RIFF Identifier
    header_raw_4b = f.read(4)
    riff_idt = header_raw_4b.decode('ascii')
    print('RIFF Identifyier:', riff_idt)

    # read chank size
    header_raw_4b = f.read(4)
    #print(header_raw_4b)
    chank_size = struct.unpack("<I", header_raw_4b)[0]
    print('chank size:', chank_size, 'Bytes')

    # read format
    header_raw_4b = f.read(4)
    file_format = header_raw_4b.decode('ascii')
    print('file format:', file_format)

    # fmt identifier
    header_raw_4b = f.read(4)
    fmt_idt = header_raw_4b.decode('ascii')
    print('fmt identifier:', fmt_idt)

    # bytes of fmt chank
    header_raw_4b = f.read(4)
    byte_size_of_chank = struct.unpack("<I", header_raw_4b)[0]
    print('byte_size_of_chank:', byte_size_of_chank)

    # format
    header_raw_2b = f.read(2)
    voice_format = struct.unpack("<H", header_raw_2b)[0]
    print('voice_format:', voice_format)

    # n of channels
    header_raw_2b = f.read(2)
    n_of_channels = struct.unpack("<H", header_raw_2b)[0]
    print('n_of_channels:', n_of_channels, 'ch')

    # sampling freq
    header_raw_4b = f.read(4)
    sampling_ferq = struct.unpack("<I", header_raw_4b)[0]
    print('sampling_ferq:', sampling_ferq, 'Hz')

    # Byte per sec
    header_raw_4b = f.read(4)
    voice_bps = struct.unpack("<I", header_raw_4b)[0]
    print('voice_Bps:', voice_bps, 'Bps')

    # block size
    header_raw_2b = f.read(2)
    block_size = struct.unpack("<H", header_raw_2b)[0]
    print('block_size:', block_size)

    # bits per sample
    header_raw_2b = f.read(2)
    bits_per_sample = struct.unpack("<H", header_raw_2b)[0]
    print('bits_per_sample:', bits_per_sample)

    # skip extention params(2B)
    # skip extention params(NB)

    # sub chank identifier 
    header_raw_4b = f.read(4)
    chank_idt = header_raw_4b.decode('ascii')
    print('sub chank identifier:', chank_idt)

    # size of sub chank 
    header_raw_4b = f.read(4)

    print('chank identifier:', chank_idt)
    size_of_sub_chank = struct.unpack("<I", header_raw_4b)[0]
    print('sub chank size:', size_of_sub_chank, 'Bytes')
    return size_of_sub_chank

def get_n_of_samples(target_freq):
    n = int(SAMPLING_FREQUENCY / target_freq)
    return n

def read_one_sample(fp):
    # read 1 sample flame
    sample_raw_4b = fp.read(4)
    sample_l, sample_r = struct.unpack("<hh", sample_raw_4b)
    monoral = int((sample_l + sample_r) / 2)
    return(monoral)

def read_samples(fp, size):
   buffer=[]
   for i in range(size):
        sample = read_one_sample(fp)
        buffer.append(sample)
   return buffer


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)

def dft(freq, sampling_freq, buffer):
    freq_rel = 0
    freq_im = 0
    for i in range(len(buffer)):
        signal = buffer[i]
        v_sin, v_cos = get_sin_cos(freq, sampling_freq, i)
        freq_rel += signal * v_cos
        freq_im += signal * v_sin
    spectrum = math.sqrt(freq_rel**2 + freq_im**2) / WINDOW_SIZE
    return spectrum

OVERLAP_RATE = 50  # 50%

with open(IN_FILE,'rb') as fp_r, open(OUT_FILE,'w') as fp_w:
    size_of_sub_chank = read_WAV_headers(fp_r)   
    n_of_samples = int(size_of_sub_chank / 2 / 2)   # 2B/ch  2CH
    
    overlap_size = int(WINDOW_SIZE *  (100 - OVERLAP_RATE) / 100)
    n_of_blocks = int((n_of_samples - WINDOW_SIZE) / overlap_size)
    for freq in range(MIN_ANALYZE_FREQ, MAX_ANALYZE_FREQ+1, FREQUENCY_RESOLUTION):
        fp_w.write(f'{freq},')
    fp_w.write('\n')

    for nth_block in range(n_of_blocks):
        if nth_block == 0:
            buffer = read_samples(fp_r, WINDOW_SIZE)
        else:
            buffer = buffer[:overlap_size]  # drop prev data and keep overlap
            new_buffer = read_samples(fp_r, WINDOW_SIZE - overlap_size)
            buffer += new_buffer
            print(f'{nth_block}')
            for freq in range(MIN_ANALYZE_FREQ, MAX_ANALYZE_FREQ+1, FREQUENCY_RESOLUTION):
                spectrum = dft(freq, SAMPLING_FREQUENCY, buffer)
                fp_w.write(f'{spectrum},')
            fp_w.write('\n')