ついにWATLABブログから書籍「いきなりプログラミングPython」が発売しました！
スペクトログラムは時間×周波数×振幅と数多くの情報を得ることが可能な便利な分析手法です。ここでは信号処理の技術でスペクトログラムのピークを自動検出するPythonコードを紹介します。
こんにちは。wat(@watlablog)です。ここではスペクトログラム（カラーマップ）から任意個数のピークを抽出する方法を紹介します！
スペクトログラムの概要とピーク検出方針
スペクトログラムとは？
本記事ではスペクトログラム(Spectrogram)という信号処理分析手法を扱います。
スペクトログラムとは、横軸に時間、縦軸に周波数、奥行き（色成分）に振幅レベルをとったプロットの事で、一度に信号の時間変化、周波数変化、振幅レベル変化を読み取ることができ、特に振動騒音分野で重宝されています（以下に参考図を示します）。
スペクトログラムの概要や、Pythonでスペクトログラムをプロットする方法の詳細は「Pythonで音のSTFT計算を自作!スペクトログラム表示する方法」に記載していますので、是非ご覧下さい。
ピーク自動検出のメリット
周波数分析をするとある周波数で卓越するピーク（峰と呼ぶ人も）が出てくる事が多く、このピークの変化を追う事で計測された情報から様々な現象を分析する事が可能です。
当WATLABブログでは「PythonでFFT波形から任意個数のピークを自動検出する方法」で単純な周波数波形からピークを自動検出する方法をPythonコードで書いてみました。
基本的なメリットはこちらの記事で紹介している内容と同じで、ピークを手動で読み取るよりは自動でリスト化してくれた方が生産性が向上するというものです。
しかしながらスペクトログラムは周波数波形の時間変化を表しているため、手動でピークを読み取ろうとするとかなり大変な作業になるはずです。
そこで、本記事では周波数波形だけでなく、スペクトログラムにおいてもピーク自動検出機能を実装してみようと思い書き始めました。
スペクトログラムのピーク検出方針(大きい順に抽出する方法）
ピークを検出する目的は様々ですが、ここではスペクトログラムから大きい順にピークを抽出する方法を考えていきます。
他にも、N番目のピークのみを抽出し、回転パルス無しにトラッキング解析をしたり…とか、ピーク活用は目的に応じて色々ありますね。
既存のピーク検出アルゴリズムを使う
スペクトログラムからピークを検出する…といっても基本的には先に紹介した記事にあるように、既存のピーク検出アルゴリズムをそのまま流用します。
既存のピーク検出アルゴリズム（ここではscipyのargrelmax等）はorderの指定が可能でノイズ対策がとれていたり、任意個数を指定できたり、自分で微分計算から作るとやらなければならない手間を大幅に削減できるというのが最大の理由です。
スペクトログラムは以下の図のように縦方向に時間成分でスライスすれば単一の周波数波形になるので、この形にすることができれば既存のピーク検出アルゴリズムを適用することが可能です。
ピーク情報をストックする
単一の周波数波形のみを検出しただけではスペクトログラムのピーク検出をしたとは言えません。
そのため、下図に示すように初期時間\(t_{1}\)から順番に\(t_{2}\)…そして最終時刻\(t_{m}\)までピーク検出を行います。
検出したピークは振幅レベル(Peak), 周波数(Freq), 時間(t)の情報としてそれぞれ2次元配列に格納していきます（下図を参考のこと）。
ピーク情報を平坦化する
上記方法でピーク情報をストックしていくと、非常に大量のピークが検出されてしまうことがあります。ピーク検出アルゴリズムは重要なピークのみを抽出する事が必要です。
ここでは振幅レベルでソートできる形にするため、先ほどの2次元配列を平坦化（1次元配列）にする方針をとります（以下参考図）。
ソートして任意個数を抽出する
ここまでできれば、あとは振幅レベル(Peak)でソートし、周波数(Freq)、時間(t)情報も並べ替えした指標で抽出してくれば、それぞれの配列の先頭から任意個数を抽出するのは簡単になります。
スペクトログラムからピーク検出するPythonコード
ピーク検出の関数
今回メインとなるスペクトログラムからピークを検出するPythonコードをdef関数として以下に示します。考案したアルゴリズムは上記説明の通りですので、コード内コメントを参照して頂ければおそらく読み解けると思います。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 |
import numpy as np from scipy import signal # 波形(x, y)からn個のピークを幅wで検出する関数(xは0から始まる仕様） def findpeaks(dx, y, n, w): index_all = list(signal.argrelmax(y, order=w)) # scipyのピーク検出 index = [] # ピーク指標の空リスト peaks = [] # ピーク値の空リスト # n個分のピーク情報(指標、値）を格納 for i in range(n): # n個のピークに満たない場合は途中でループを抜ける（エラー処理） if i >= len(index_all[0]): break index.append(index_all[0][i]) peaks.append(y[index_all[0][i]]) # 個数の足りない分を0で埋める（エラー処理） if len(index) != n: index = index + ([0] * (n - len(index))) peaks = peaks + ([0] * (n - len(peaks))) index = np.array(index) * dx # xの分解能dxをかけて指標を物理軸に変換 peaks = np.array(peaks) return index, peaks # スペクトログラムからピークを検出する関数 # fft_array=スペクトログラム（転置前） # dt, df=スペクトログラムの時間分解能, 周波数分解能 # num_peaks=1つの周波数軸で検出するピーク数 # w=ノイズ対策用の幅（order） # max_peaks=最終的にスペクトログラムから抽出するピーク数（振幅の大きい順） def findpeaks_2d(fft_array, dt, df, num_peaks, w, max_peaks): # ピーク情報を初期化する time_index = np.zeros((len(fft_array), num_peaks)) freq_index = np.zeros((len(fft_array), num_peaks)) freq_peaks = np.zeros((len(fft_array), num_peaks)) # 各周波数軸毎にピークを検出する for i in range(len(fft_array)): index, peaks = findpeaks(df, fft_array[i], n=num_peaks, w=w) # ピーク検出 freq_peaks[i] = peaks # 検出したピーク値(振幅)を格納 freq_index[i] = index # 検出したピーク位置(周波数)を格納 time_index[i] = np.full(num_peaks, i) * dt # 検出したピーク位置(時間)を格納 # 平坦化する freq_peaks = freq_peaks.ravel() freq_index = freq_index.ravel() time_index = time_index.ravel() # ピークの大きい順（降順）にソートする freq_peaks_sort = np.sort(freq_peaks)[::-1] freq_index_sort = freq_index[np.argsort(freq_peaks)[::-1]] time_index_sort = time_index[np.argsort(freq_peaks)[::-1]] return freq_index_sort[:max_peaks], freq_peaks_sort[:max_peaks], time_index_sort[:max_peaks] |
全コード（コピペ用：単一ピークを検出）
以下にコピペしてすぐに動かせるように全コードを示します。.pyファイルのあるディレクトリに「sample.wav」（ここでは5[s]のデータ）というwavファイルを置くことでそのまま動きます。
その他の関数群は「Pythonで音のSTFT計算を自作!スペクトログラム表示する方法」で使っていたものと同じですので、もし詳細で不明な所があればこの記事を参照下さい。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 |
import numpy as np from scipy import signal from scipy import fftpack import soundfile as sf from matplotlib import pyplot as plt # オーバーラップをかける関数 def ov(data, samplerate, Fs, overlap): Ts = len(data) / samplerate # 全データ長 Fc = Fs / samplerate # フレーム周期 x_ol = Fs * (1 - (overlap / 100)) # オーバーラップ時のフレームずらし幅 N_ave = int((Ts - (Fc * (overlap / 100))) / (Fc * (1 - (overlap / 100)))) # 抽出するフレーム数（平均化に使うデータ個数） array = [] # 抽出したデータを入れる空配列の定義 # forループでデータを抽出 for i in range(N_ave): ps = int(x_ol * i) # 切り出し位置をループ毎に更新 array.append(data[ps:ps + Fs:1]) # 切り出し位置psからフレームサイズ分抽出して配列に追加 final_time = (ps + Fs)/samplerate # 切り出したデータの最終時刻 return array, N_ave, final_time # オーバーラップ抽出されたデータ配列とデータ個数、最終時間を戻り値にする # ハニング窓をかける関数（振幅補正係数計算付き） def hanning(data_array, Fs, N_ave): han = signal.hann(Fs) # ハニング窓作成 acf = 1 / (sum(han) / Fs) # 振幅補正係数(Amplitude Correction Factor) # オーバーラップされた複数時間波形全てに窓関数をかける for i in range(N_ave): data_array[i] = data_array[i] * han # 窓関数をかける return data_array, acf # dB(デシベル）演算する関数 def db(x, dBref): y = 20 * np.log10(x / dBref) # 変換式 return y # dB値を返す # 聴感補正(A補正)する関数 def aweightings(f): if f[0] == 0: f[0] = 1 ra = (np.power(12194, 2) * np.power(f, 4))/\ ((np.power(f, 2) + np.power(20.6, 2)) * np.sqrt((np.power(f, 2) + np.power(107.7, 2)) * (np.power(f, 2) + np.power(737.9, 2))) * (np.power(f, 2) + np.power(12194, 2))) a = 20 * np.log10(ra) + 2.00 return a # 平均化FFTする関数 def fft_ave(data_array, samplerate, Fs, N_ave, acf): fft_array = [] fft_axis = np.linspace(0, samplerate, Fs) # 周波数軸を作成 a_scale = aweightings(fft_axis) # 聴感補正曲線を計算 # FFTをして配列にdBで追加、窓関数補正値をかけ、(Fs/2)の正規化を実施。 for i in range(N_ave): fft_array.append(db (acf * np.abs(fftpack.fft(data_array[i]) / (Fs / 2)) , 2e-5)) fft_array = np.array(fft_array) + a_scale # 型をndarrayに変換しA特性をかける fft_mean = np.mean(fft_array, axis=0) # 全てのFFT波形の平均を計算 return fft_array, fft_mean, fft_axis # 波形(x, y)からn個のピークを幅wで検出する関数(xは0から始まる仕様） def findpeaks(dx, y, n, w): index_all = list(signal.argrelmax(y, order=w)) # scipyのピーク検出 index = [] # ピーク指標の空リスト peaks = [] # ピーク値の空リスト # n個分のピーク情報(指標、値）を格納 for i in range(n): # n個のピークに満たない場合は途中でループを抜ける（エラー処理） if i >= len(index_all[0]): break index.append(index_all[0][i]) peaks.append(y[index_all[0][i]]) # 個数の足りない分を0で埋める（エラー処理）-->最小値が0以下になるようなデータの場合は検討要 if len(index) != n: index = index + ([0] * (n - len(index))) peaks = peaks + ([0] * (n - len(peaks))) index = np.array(index) * dx # xの分解能dxをかけて指標を物理軸に変換 peaks = np.array(peaks) return index, peaks # スペクトログラムからピークを検出する関数 # fft_array=スペクトログラム（転置前） # dt, df=スペクトログラムの時間分解能, 周波数分解能 # num_peaks=1つの周波数軸で検出するピーク数 # w=ノイズ対策用の幅（order） # max_peaks=最終的にスペクトログラムから抽出するピーク数（振幅の大きい順） def findpeaks_2d(fft_array, dt, df, num_peaks, w, max_peaks): # ピーク情報を初期化する time_index = np.zeros((len(fft_array), num_peaks)) freq_index = np.zeros((len(fft_array), num_peaks)) freq_peaks = np.zeros((len(fft_array), num_peaks)) # 各周波数軸毎にピークを検出する for i in range(len(fft_array)): index, peaks = findpeaks(df, fft_array[i], n=num_peaks, w=w) # ピーク検出 freq_peaks[i] = peaks # 検出したピーク値(振幅)を格納 freq_index[i] = index # 検出したピーク位置(周波数)を格納 time_index[i] = np.full(num_peaks, i) * dt # 検出したピーク位置(時間)を格納 # 平坦化する freq_peaks = freq_peaks.ravel() freq_index = freq_index.ravel() time_index = time_index.ravel() # ピークの大きい順（降順）にソートする freq_peaks_sort = np.sort(freq_peaks)[::-1] freq_index_sort = freq_index[np.argsort(freq_peaks)[::-1]] time_index_sort = time_index[np.argsort(freq_peaks)[::-1]] return freq_index_sort[:max_peaks], freq_peaks_sort[:max_peaks], time_index_sort[:max_peaks] path = 'sample.wav' #ファイルパスを指定 data, samplerate = sf.read(path) #wavファイルを読み込む x = np.arange(0, len(data)) / samplerate #波形生成のための時間軸の作成 # Fsとoverlapでスペクトログラムの分解能を調整する。 Fs = 4096 # フレームサイズ overlap = 90 # オーバーラップ率 # オーバーラップ抽出された時間波形配列 time_array, N_ave, final_time = ov(data, samplerate, Fs, overlap) # ハニング窓関数をかける time_array, acf = hanning(time_array, Fs, N_ave) # FFTをかける fft_array, fft_mean, fft_axis = fft_ave(time_array, samplerate, Fs, N_ave, acf) # ピーク検出する spec_dt = final_time / N_ave freq_index, freq_peaks, time_index = findpeaks_2d(fft_array, spec_dt, fft_axis[1], num_peaks=100, w=10, max_peaks=1) print(freq_peaks) print(freq_index) print(time_index) # スペクトログラムで縦軸周波数、横軸時間にするためにデータを転置 fft_array = fft_array.T # ここからグラフ描画 # グラフをオブジェクト指向で作成する。 fig = plt.figure() ax1 = fig.add_subplot(111) # データをプロットする。 im = ax1.imshow(fft_array, vmin=0, vmax=80, extent=[0, final_time, 0, samplerate], aspect='auto', cmap='jet') # カラーバーを設定する。 cbar = fig.colorbar(im) cbar.set_label('SPL [dBA]') # ピークを描画する。 ax1.scatter(time_index, freq_index, s=250, facecolor='None', linewidth=3, edgecolors='yellow') for i in range(len(time_index)): plt.text(time_index[i], freq_index[i] + 150, str(round(freq_peaks[i], 2)), fontsize=10, color='black', backgroundcolor='white', horizontalalignment='center') # 軸設定する。 ax1.set_xlabel('Time [s]') ax1.set_ylabel('Frequency [Hz]') # スケールの設定をする。 ax1.set_xticks(np.arange(0, 50, 1)) ax1.set_yticks(np.arange(0, 20000, 500)) ax1.set_xlim(0, 5) ax1.set_ylim(0, 2000) # グラフを表示する。 plt.show() plt.close() |
上記コードを実行すると、スペクトログラム内で最も最大振幅となる点を教えてくれます。
グラフとしては以下のようなプロットがでます。
コンソールには以下のようにピーク情報が表示されます。
|
1 2 3 |
[77.41831688] # 振幅レベル[dB] [1023.07692308] # 周波数[Hz] [2.12541551] # 時間[s] |
複数ピークを検出する方法
先ほどは単一のピークでしたが、max_peaksの値を変更することで複数ピークをスペクトログラム全体から振幅レベルの大きい順に検出することが出来ます。
|
1 |
freq_index, freq_peaks, time_index = findpeaks_2d(fft_array, spec_dt, fft_axis[1], num_peaks=100, w=10, max_peaks=10) |
以下はmax_peaks=10, 50の比較です。数値を増やしていけば検出するピークも増えていきます。
コンソールには検出したピーク情報が配列で表示されます。
|
1 2 3 4 5 6 7 8 |
[77.41831688 77.15518915 76.99655167 76.23593015 76.22668513 75.94165624 75.76663923 75.38825859 75.22817809 75.16880769] [1023.07692308 1023.07692308 1023.07692308 1023.07692308 1130.76923077 1023.07692308 1130.76923077 1130.76923077 990.76923077 990.76923077] [2.12541551 2.11596922 2.1348618 2.10652292 2.6638541 2.14430809 2.65440781 2.67330039 1.8136879 1.82313419] |
録音データからピーク検出するコード
スペクトログラムは録音データに対して計算する需要も多くあるようです。そのため以下にPyAudioを使った録音版のスペクトログラムピーク検出コードを載せておきます。
PyAudioによる録音については「PythonのPyAudioで音声録音をする簡単な方法」という記事に詳細を載せていますので、是非ご覧下さい。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 |
import numpy as np from scipy import signal from scipy import fftpack from matplotlib import pyplot as plt import pyaudio def record(index, samplerate, fs, time): pa = pyaudio.PyAudio() data = [] dt = 1 / samplerate # ストリームの開始 stream = pa.open(format=pyaudio.paInt16, channels=1, rate=samplerate, input=True, input_device_index=index, frames_per_buffer=fs) # フレームサイズ毎に音声を録音していくループ for i in range(int(((time / dt) / fs))): frame = stream.read(fs) data.append(frame) # ストリームの終了 stream.stop_stream() stream.close() pa.terminate() # データをまとめる処理 data = b"".join(data) # データをNumpy配列に変換 data = np.frombuffer(data, dtype="int16") / float((np.power(2, 16) / 2) - 1) return data, i # オーバーラップをかける関数 def ov(data, samplerate, Fs, overlap): Ts = len(data) / samplerate # 全データ長 Fc = Fs / samplerate # フレーム周期 x_ol = Fs * (1 - (overlap / 100)) # オーバーラップ時のフレームずらし幅 N_ave = int((Ts - (Fc * (overlap / 100))) / (Fc * (1 - (overlap / 100)))) # 抽出するフレーム数（平均化に使うデータ個数） array = [] # 抽出したデータを入れる空配列の定義 # forループでデータを抽出 for i in range(N_ave): ps = int(x_ol * i) # 切り出し位置をループ毎に更新 array.append(data[ps:ps + Fs:1]) # 切り出し位置psからフレームサイズ分抽出して配列に追加 final_time = (ps + Fs)/samplerate # 切り出したデータの最終時刻 return array, N_ave, final_time # オーバーラップ抽出されたデータ配列とデータ個数、最終時間を戻り値にする # ハニング窓をかける関数（振幅補正係数計算付き） def hanning(data_array, Fs, N_ave): han = signal.hann(Fs) # ハニング窓作成 acf = 1 / (sum(han) / Fs) # 振幅補正係数(Amplitude Correction Factor) # オーバーラップされた複数時間波形全てに窓関数をかける for i in range(N_ave): data_array[i] = data_array[i] * han # 窓関数をかける return data_array, acf # dB(デシベル）演算する関数 def db(x, dBref): y = 20 * np.log10(x / dBref) # 変換式 return y # dB値を返す # 聴感補正(A補正)する関数 def aweightings(f): if f[0] == 0: f[0] = 1 ra = (np.power(12194, 2) * np.power(f, 4))/\ ((np.power(f, 2) + np.power(20.6, 2)) * np.sqrt((np.power(f, 2) + np.power(107.7, 2)) * (np.power(f, 2) + np.power(737.9, 2))) * (np.power(f, 2) + np.power(12194, 2))) a = 20 * np.log10(ra) + 2.00 return a # 平均化FFTする関数 def fft_ave(data_array, samplerate, Fs, N_ave, acf): fft_array = [] fft_axis = np.linspace(0, samplerate, Fs) # 周波数軸を作成 a_scale = aweightings(fft_axis) # 聴感補正曲線を計算 # FFTをして配列にdBで追加、窓関数補正値をかけ、(Fs/2)の正規化を実施。 for i in range(N_ave): fft_array.append(db (acf * np.abs(fftpack.fft(data_array[i]) / (Fs / 2)) , 2e-5)) fft_array = np.array(fft_array) + a_scale # 型をndarrayに変換しA特性をかける fft_mean = np.mean(fft_array, axis=0) # 全てのFFT波形の平均を計算 return fft_array, fft_mean, fft_axis # 波形(x, y)からn個のピークを幅wで検出する関数(xは0から始まる仕様） def findpeaks(dx, y, n, w): index_all = list(signal.argrelmax(y, order=w)) # scipyのピーク検出 index = [] # ピーク指標の空リスト peaks = [] # ピーク値の空リスト # n個分のピーク情報(指標、値）を格納 for i in range(n): # n個のピークに満たない場合は途中でループを抜ける（エラー処理） if i >= len(index_all[0]): break index.append(index_all[0][i]) peaks.append(y[index_all[0][i]]) # 個数の足りない分を0で埋める（エラー処理）-->最小値が0以下になるようなデータの場合は検討要 if len(index) != n: index = index + ([0] * (n - len(index))) peaks = peaks + ([0] * (n - len(peaks))) index = np.array(index) * dx # xの分解能dxをかけて指標を物理軸に変換 peaks = np.array(peaks) return index, peaks # スペクトログラムからピークを検出する関数 # fft_array=スペクトログラム（転置前） # dt, df=スペクトログラムの時間分解能, 周波数分解能 # num_peaks=1つの周波数軸で検出するピーク数 # w=ノイズ対策用の幅（order） # max_peaks=最終的にスペクトログラムから抽出するピーク数（振幅の大きい順） def findpeaks_2d(fft_array, dt, df, num_peaks, w, max_peaks): # ピーク情報を初期化する time_index = np.zeros((len(fft_array), num_peaks)) freq_index = np.zeros((len(fft_array), num_peaks)) freq_peaks = np.zeros((len(fft_array), num_peaks)) # 各周波数軸毎にピークを検出する for i in range(len(fft_array)): index, peaks = findpeaks(df, fft_array[i], n=num_peaks, w=w) # ピーク検出 freq_peaks[i] = peaks # 検出したピーク値(振幅)を格納 freq_index[i] = index # 検出したピーク位置(周波数)を格納 time_index[i] = np.full(num_peaks, i) * dt # 検出したピーク位置(時間)を格納 # 平坦化する freq_peaks = freq_peaks.ravel() freq_index = freq_index.ravel() time_index = time_index.ravel() # ピークの大きい順（降順）にソートする freq_peaks_sort = np.sort(freq_peaks)[::-1] freq_index_sort = freq_index[np.argsort(freq_peaks)[::-1]] time_index_sort = time_index[np.argsort(freq_peaks)[::-1]] return freq_index_sort[:max_peaks], freq_peaks_sort[:max_peaks], time_index_sort[:max_peaks] # Fsとoverlapでスペクトログラムの分解能を調整する。 Fs = 4096 # フレームサイズ overlap = 90 # オーバーラップ率 # 録音して波形を作る（indexはPCによる） index = 0 samplerate = 44100 t_max = 5 data, i = record(index, samplerate, Fs, t_max) x = np.arange(0, Fs * (i+1) * (1 / samplerate), 1 / samplerate) # オーバーラップ抽出された時間波形配列 time_array, N_ave, final_time = ov(data, samplerate, Fs, overlap) # ハニング窓関数をかける time_array, acf = hanning(time_array, Fs, N_ave) # FFTをかける fft_array, fft_mean, fft_axis = fft_ave(time_array, samplerate, Fs, N_ave, acf) # ピーク検出する spec_dt = final_time / N_ave freq_index, freq_peaks, time_index = findpeaks_2d(fft_array, spec_dt, fft_axis[1], num_peaks=100, w=10, max_peaks=1) print(freq_peaks) print(freq_index) print(time_index) # スペクトログラムで縦軸周波数、横軸時間にするためにデータを転置 fft_array = fft_array.T # ここからグラフ描画 # グラフをオブジェクト指向で作成する。 fig = plt.figure() ax1 = fig.add_subplot(111) # データをプロットする。 im = ax1.imshow(fft_array, vmin=0, vmax=80, extent=[0, final_time, 0, samplerate], aspect='auto', cmap='jet') # カラーバーを設定する。 cbar = fig.colorbar(im) cbar.set_label('SPL [dBA]') # ピークを描画する。 ax1.scatter(time_index, freq_index, s=250, facecolor='None', linewidth=3, edgecolors='yellow') for i in range(len(time_index)): plt.text(time_index[i], freq_index[i] + 150, str(round(freq_peaks[i], 2)), fontsize=10, color='black', backgroundcolor='white', horizontalalignment='center') # 軸設定する。 ax1.set_xlabel('Time [s]') ax1.set_ylabel('Frequency [Hz]') # スケールの設定をする。 ax1.set_xticks(np.arange(0, 50, 1)) ax1.set_yticks(np.arange(0, 20000, 500)) ax1.set_xlim(0, 5) ax1.set_ylim(0, 2000) # グラフを表示する。 plt.show() plt.close() |
以下が例です。口笛を吹いてみました。
まとめ
本記事では過去に当WATLABブログで紹介したスペクトログラムの概要を再度説明し、スペクトログラムからピークを自動検出することの有用性を紹介しました。
記事の中盤では考案したピーク検出アルゴリズムを図解で説明しました。
最後に考案したアルゴリズムをPythonコード化し、サンプルのwavファイルを使って実際にピーク検出を行い、パラメータ違いで任意個数のピークを検出できることを確認しました。
本ピーク検出アルゴリズムは一例に過ぎません。是非各自の目的に合わせてカスタマイズしてみて下さい（特に当ブログのコードをコピペして使ったり2次利用するのに使用許諾とかは必要ありません）。
スペクトログラムという2次元データから任意個数のピークを抽出することが出来ました！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！
いつも参考にさせていただいております。wat様のscipy.signal.argrelmaxを使ってFFTのピーク値を検出するという試みのサイトも拝見させていただきまして、どちらにコメントさせていただこうかと迷いましたがこちらで失礼いたします。
私もピーク値を検出するプログラムを作成しているのですが、なかなかうまくいきません。scipyのsignalのなかにはいっているspectrogram()を用いてスペクトログラムを作成し、最終的にはwat様の本ページのようなグラフ内ピーク値の描画ができないかと検討しています。
まだPythonを初めて間もなく、あまりにも抽象的で初歩的な質問で申し訳ないのですが、そもそもspectrogram()を用いたうえでピーク値を特定することは可能でしょうか。
ご教授いただければと思います。よろしくお願いいたします。
いつもご訪問ありがとうございます。
scipyのspectrogramは使った事ないのですが、試しに
f,t,Sxx = signal.spectrogram(data, Fs, nperseg=512)
print(Sxx)
と書いてみますと、Sxxは2次元の結果が返ってきます。
公式リファレンス
signal.spectrogram()
を確認すると、fとtはそれぞれ周波数軸と時間軸のようです。
本ページに記載のfindpeaks_2d()関数もスペクトログラムを2次元配列の引数としてデータを渡しているので、本質的には同じデータの形をしていますので、signal.spectrogram()の結果からもピーク検出はできる事になります。
mode=magnitudeで振幅成分を出力する事もできそうですが、おそらくこの値はdB変換等はしていない値だと思いますので、
綺麗にグラフ上で次数を表示させるには、色レベルをログスケールにするか、dB変換を行うのが良いと考えます。
ただ、少し公式リファレンスを良く読んでsignal.spectrogram()の動作を良く理解してから使わないと目的の値を得るのは難しいと思います。
まずはnumpy等で理論波形を作り、意図したスペクトログラムが表示されるプログラムを作成した後で、ピーク検出を検討するのが良いと思いますが、
現在なめこ様は既にスペクトログラム表示まではできていますでしょうか？
返信ありがとうございます。
スペクトログラム表示まではできております。
先月から始めたPythonで知識もおぼつかず、見当違いなことをしてしまっているのだと思います。
スペクトログラムを行う信号はpyaudioを使った録音で波形をとる形をとっております。そのうえでargrelmax()を使ってピーク値を出そうと試みたのですがうまくいかず…という感じで昨日が終わりました。
既にTwitterでのやりとりで解決しましたが、他の方のためにこちらにも記録しておきます。
■scipy.signal.spectrogramについて
以下の記事にトライ結果を追記しました。
「Pythonで音のSTFT計算を自作!スペクトログラム表示する方法」
→STFTからスペクトログラム表示までは問題なくできましたが、振幅成分に不明点あり使用には注意が必要です（記事内の自作コードであればそれらの不明点はありません）。
■ピーク検出の方法とエラー対処について
エラーの原因はピーク検出時のインデックスに不備があったからでした。
「PythonでFFT波形から任意個数のピークを自動検出する方法」の記事に基本のピーク検出関数にエラー処理を入れる事で対応可能でした。
最終的に録音データのスペクトログラムからピーク検出するコードはこちらにまとめました。
「Pythonでスペクトログラムからピーク値を任意数抽出する方法」
なめこ様のおかげでプログラムのバグ出しができ、記事を改善する事ができました。
ありがとうございます。
今後もよろしくお願い致します。