音の時間波形の包絡線を抽出や振幅変調(AM)の復調などで利用されるヒルベルト変換をPythonで実装し、いくつかの簡単な信号に対して適用してみました。
ヒルベルト変換の概要
ヒルベルト変換の仕組みや実装方法については、以下のWebページがわかりやすいです。
小野測器 基礎からの周波数分析(29)-「ヒルベルト変換と解析信号」
この記事に従い,
- フーリエ変換を用いて解析信号を算出
- 解析信号から包絡線,瞬時位相,瞬時周波数を抽出
する処理を実装します.
Python実装
入力信号→解析信号の算出
def hilbert_transform(x): n = x.shape[0] x_f = np.fft.fft(x, n) x_f[1:n//2+1] = 2*x_f[1:n//2+1] x_f[n//2+1:] = 0.0 z = np.fft.ifft(x_f, n) return z
解析信号→包絡線,瞬時位相,瞬時周波数の抽出
def compute_instantaneous_params(z, sr): x, x_h = z.real, z.imag envelope = np.sqrt(x**2 + x_h**2) phase = np.angle(z) #arctan2(x_h, x) freq = np.diff(np.unwrap(phase))/(2*np.pi)*sr return envelope, phase, freq
信号への適用例
正弦波への適用例
基本周波数が 200Hz の正弦波に適用した結果です。
解析信号の虚部が、元信号から90度位相がずれた信号になっていることが確認できます。
単純な正弦波では、包絡線は一定の値となっています。
[px]
振幅変調された信号
次に、先ほどの正弦波に対して,10Hz の正弦波による振幅変調*1を適用した信号を分析してみます。
結果として、10Hz の正弦波が包絡線として抽出できていることを確認しました。
[px]
減衰する正弦波
[px]
周波数が時間変化する正弦波(チャープ信号)
瞬時周波数が上昇していることが確認できました。
[px]
プロット用Python スクリプト
import numpy as np import matplotlib.pyplot as plt # https://www.onosokki.co.jp/HP-WK/eMM_back/emm180.pdf def hilbert_transform(x): n = x.shape[0] x_f = np.fft.fft(x, n) x_f[1:n//2+1] = 2*x_f[1:n//2+1] x_f[n//2+1:] = 0.0 z = np.fft.ifft(x_f, n) return z def compute_instantaneous_params(z, sr): x, x_h = z.real, z.imag envelope = np.sqrt(x**2 + x_h**2) phase = np.angle(z) #arctan2(x_h, x) freq = np.diff(np.unwrap(phase))/(2*np.pi)*sr return envelope, phase, freq sr = 48000 # Sampling rate duration = 0.3 # [sec] fo = 200 fam = 10 t = np.arange(0, sr*duration) / sr ## 正弦波 x = np.sin(2*np.pi*fo*t) ## 振幅変調 #xam = 0.5*np.sin(2*np.pi*fam*t) + 0.5 #x = xam * x # ## 減衰する正弦波 #x *= np.exp(-10*t) ## 周波数が時間変化する正弦波 #from scipy.signal import chirp #x = chirp(t, f0=100, f1=1000, t1=duration, method='log') z = hilbert_transform(x) envelope, instant_phase, instant_freq = compute_instantaneous_params(z, sr) plt.subplot(4,1,1) plt.plot(x, label="x") plt.legend(loc="lower right") plt.grid() plt.xlim([0,sr*duration]) plt.ylim([-1.5,1.5]) plt.subplot(4,1,2) plt.plot(z.real, label="z.real", alpha=0.5) plt.plot(z.imag, label="z.imag", alpha=0.5) plt.plot(envelope, label="envelope") plt.legend(loc="lower right") plt.grid() plt.xlim([0,sr*duration]) plt.ylim([-1.5,1.5]) plt.subplot(4,1,3) plt.plot(instant_phase, label="instant. phase") plt.legend(loc="lower right") plt.grid() plt.xlim([0,sr*duration]) plt.ylim([-np.pi,np.pi]) plt.subplot(4,1,4) plt.plot(instant_freq, label="instant. freq") plt.legend(loc="lower right") plt.xlim([0,sr*duration]) plt.ylim([100,1000]) plt.grid() plt.tight_layout() plt.show()
*1:平均0.5, 最大値1.0, 最小値0.0