Wizard Notes

Python, JavaScript を使った音楽信号分析の技術録、作曲活動に関する雑記

Python:ヒルベルト変換で様々な信号の包絡線,瞬時位相,瞬時周波数を抽出

音の時間波形の包絡線を抽出や振幅変調(AM)の復調などで利用されるヒルベルト変換Pythonで実装し、いくつかの簡単な信号に対して適用してみました。

ヒルベルト変換の概要

ヒルベルト変換の仕組みや実装方法については、以下のWebページがわかりやすいです。

小野測器 基礎からの周波数分析(29)-「ヒルベルト変換と解析信号」

この記事に従い,

  1. フーリエ変換を用いて解析信号を算出
  2. 解析信号から包絡線,瞬時位相,瞬時周波数を抽出

する処理を実装します.

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