python スペクトルピーク検出

前提・実現したいこと

録音した音声のスペクトルピークを検出したい

発生している問題・エラーメッセージ

Traceback (most recent call last):
  File "/Users/matsumura/PycharmProjects/wave/FFT.py", line 114, in <module>
    aaa, bbb = find_peaks(xs, 256)
  File "/Users/matsumura/PycharmProjects/wave/FFT.py", line 49, in find_peaks
    return (result_idxs, a[result_idxs])
IndexError: arrays used as indices must be of integer (or boolean) type

該当のソースコード

python
1#!/usr/bin/env python
2# vim:fileencoding=utf-8
3
4import pyaudio as pa
5import time
6import numpy as np
7import soundfile as sf
8import sys
9import scipy.signal as sg
10import matplotlib.pyplot as plt
11
12xs = np.array([])
13
14
15def find_peaks(a, amp_thre, local_width=1, min_peak_distance=1):
16    """
17    閾値と極大・極小を判定する窓幅、ピーク間最小距離を与えて配列からピークを検出する。
18    内部的にはピーク間距離は正負で区別して算出されるため、近接した正負のピークは検出される。
19    :rtype (int, float)
20    :return tuple (ndarray of peak indices, ndarray of peak value)
21    """
22    # generate candidate indices to limit by threthold
23    idxs = np.where(np.abs(a) > amp_thre)[0]
24
25    # extend array to decide local maxima/minimum
26    idxs_with_offset = idxs + local_width
27    a_extend = np.r_[[a[0]] * local_width, a, [a[-1]] * local_width]
28
29    last_pos_peak_idx = 0
30    last_neg_peak_idx = 0
31    result_idxs = []
32
33    for i in idxs_with_offset:
34        is_local_maximum = (a_extend[i] >= 0 and
35                            a_extend[i] >= np.max(a_extend[i - local_width: i + local_width + 1]))
36        is_local_minimum = (a_extend[i] < 0 and
37                            a_extend[i] <= np.min(a_extend[i - local_width: i + local_width + 1]))
38        if (is_local_maximum or is_local_minimum):
39            if is_local_minimum:
40                if i - last_pos_peak_idx > min_peak_distance:
41                    result_idxs.append(i)
42                    last_pos_peak_idx = i
43            else:
44                if i - last_neg_peak_idx > min_peak_distance:
45                    result_idxs.append(i)
46                    last_neg_peak_idx = i
47
48    result_idxs = np.array(result_idxs) - local_width
49    return (result_idxs, a[result_idxs])
50
51
52def callback(in_data, frame_count, time_info, status):
53    global xs
54    in_float = np.frombuffer(in_data, dtype=np.int16).astype(np.float)
55    in_float[in_float > 0.0] /= float(2**15 - 1)
56    in_float[in_float <= 0.0] /= float(2**15)
57    xs = np.r_[xs, in_float]
58
59    return (in_data, pa.paContinue)
60
61
62if __name__ == "__main__":
63    # pyaudio
64    p_in = pa.PyAudio()
65    py_format = p_in.get_format_from_width(2)
66    fs = 16000
67    channels = 1
68    chunk = 1024
69    use_device_index = 0
70
71    # 入力ストリームを作成
72    in_stream = p_in.open(format=py_format,
73                          channels=channels,
74                          rate=fs,
75                          input=True,
76                          frames_per_buffer=chunk,
77                          input_device_index=use_device_index,
78                          stream_callback=callback)
79
80    in_stream.start_stream()
81
82    # input loop
83    # 何か入力したら終了
84    while in_stream.is_active():
85        c = input()
86        if c:
87            break
88        time.sleep(0.1)
89    else:
90        in_stream.stop_stream()
91        in_stream.close()
92
93    # 入力信号を保存
94    sf.write("pyaudio_output.wav", xs, fs)
95
96    p_in.terminate()
97
98
99    #デジタルフィルタを時系列信号に適用する
100    plt.close("all")
101
102    # wavファイル読み込み
103    filename = sys.argv[1]
104    wav, fs = sf.read(filename)
105
106    # # ステレオ2chの場合、LchとRchに分割
107    # wav_l = wav[:, 0]
108    # wav_r = wav[:, 1]
109
110    # 入力をモノラル化
111    # xs = (0.5 * wav_l) + (0.5 * wav_r)
112    xs = wav
113
114    aaa, bbb = find_peaks(xs, 256)
115
116    # LPF設計
117    num_tap = 1024
118    lpf_cutoff_hz = 400.0
119    lpf_cutoff = lpf_cutoff_hz / (fs/2.0)
120    win = "hann"
121    lpf = sg.firwin(num_tap, lpf_cutoff, window=win)
122
123    # 線形フィルタ適用
124    ys = sg.lfilter(lpf, [1.0], xs)
125
126    # 周波数領域フィルタ適用
127    zs = sg.fftconvolve(xs, lpf, mode="same")
128
129    # 冒頭から10秒分プロット
130    fig = plt.figure(1)
131    ax = fig.add_subplot(311)
132    ax.plot(xs[:fs*10])
133    ax.set_title("input signal")
134    ax.set_xlabel("time [pt]")
135    ax.set_ylabel("amplitude")
136
137    ax = fig.add_subplot(312)
138    ax.plot(ys[:fs*10])
139    ax.set_title("lfilter output signal")
140    ax.set_xlabel("time [pt]")
141    ax.set_ylabel("amplitude")
142
143    ax = fig.add_subplot(313)
144    ax.plot(ys[:fs*10])
145    ax.set_title("fftconvolve output signal")
146    ax.set_xlabel("time [pt]")
147    ax.set_ylabel("amplitude")
148
149    fig.set_tight_layout(True)
150
151    plt.show()