はじめに

周期性解析の手法で特許を取得して、アルゴリズムを公開して、応用先を紹介しているのですが、今回は少し変わった使い方をまとめてみたいと思います。

フーリエ変換の前処理として周期性解析を使うことによって、解析区間が倍必要になるが、信号の周波数推定の安定性を上げる方法について紹介します。

周期性について

周期性解析の手法と言えば、自己相関やフーリエ変換(以下FT)ですが、自分の特許は自己相関に近い周期性解析手法です。

下記の記事で数式の定義に触れています

下記の記事で外れ値(インパルスノイズ的な物)に強い事を紹介しています

手法

普通のFTでの周波数推定

データの切り出し
窓関数処理
FFT
指定範囲のピーク値の補完(パラボラ補完)

前処理に周期性解析の特許を使った周波数推定

データの切り出し
周期性解析処理
窓関数処理
FFT
指定範囲のピーク値の補完(パラボラ補完)

特許を使用した時間領域からの周波数推定(比較用として実験)

信号を周期解析処理にかける
指定範囲のピーク値より0.85以上の値のピーク値を取得
一番短い周期のピーク値を取得
ピーク値の補完(パラボラ補完)
周期から周波数を逆算

といった形でFFTをかける前に周期性解析処理をかけています。
解析区間が半分になるため、FFTの周波数分解能は低下します。
どちらも補完をかけて精度を向上させます。

解析データ

シミュレーションで作成します。
90Hzの波に対して、ガウスノイズとランダムインパルスノイズを加えた信号を処理します。サンプリング周波数は、20kHzで5秒間の信号としました。
サンプリング周波数に関しては、ベアリングの振動のデータに対して解析を行っていたので、その設定が残った感じですね。

結果と考察

左下のbedcmmPitchと書かれた結果が、特許を使用した時間領域からの周波数推定、
右上がF0 est - rawと書かれた結果が、普通のFTでの周波数推定、
右下のF0 est - Periodicity extractedと書かれた結果が、前処理に周期性解析の特許を使った周波数推定となっています。
特許を使用した時間領域からの周波数推定は、パラメータ調整を行い、平均値のズレ、標準偏差ともに小さくしてみました。
普通のFTでの周波数推定ですが、前処理に周期性解析の特許を使った周波数推定より、平均値のズレは小さく、標準偏差が大きくなっています。
周期性解析は周期方向への積分・平均化に近い処理となるため、推定値の分散は低下する一方、平均値に小さなオフセットが生じる結果となりました。
これは、周期性解析は平均化をしている事を示しているように感じます。

また、真ん中上の図は、生信号に窓関数をかけた信号の時間毎のFFTの結果を示した図で、真ん中下の図が、周期性解析の前処理後にFFTをした時間毎の結果となっています。
上下の図を比較すると解析の解像度が2倍違っている事がわかります。
また、上の図で見られるノイズによるまだら模様が、下の図では見えない事がわかります。細かいノイズが除去され、周期性が抽出されている事がわかります。
ランダムノイズやインパルスノイズは周期性を持たないため、周期性解析処理によって相対的に抑制されたと考えられます。
さらに、下の図では倍の周波数に少し値がある事も確認できます。これは周期性成分を抽出したことによる影響で、倍音も強調された形となっています。

特許手法のGithub

おわりに

周期性解析をピッチ検出に応用していて、やっぱりFTは強いよなと感じたので、既存のFTを強化する方法の紹介をしました。
YINに比べてインパルスノイズに強いピッチ検出方法でこちらもよかったら使ってみてください。

コード

import numpy as np
import matplotlib.pyplot as plt
from bedcmm.pattern import periodicity
import bedcmmPitch

def fft_peak_frequency(
    signal,
    fs,
    fmin,
    fmax,
    window=np.hanning
):
    N = len(signal)

    w = window(N)
    x = signal * w

    spectrum = np.fft.rfft(x)
    mag = np.abs(spectrum)

    freqs = np.fft.rfftfreq(N, 1/fs)

    idx = np.where((freqs >= fmin) & (freqs <= fmax))[0]

    local_mag = mag[idx]

    peak_local = np.argmax(local_mag)
    peak_idx = idx[peak_local]

    # parabolic interpolation
    if 1 <= peak_idx < len(mag)-1:
        alpha = mag[peak_idx - 1]
        beta  = mag[peak_idx]
        gamma = mag[peak_idx + 1]

        denom = (alpha - 2*beta + gamma)

        if abs(denom) > 1e-12:
            p = 0.5 * (alpha - gamma) / denom
        else:
            p = 0.0
    else:
        p = 0.0

    refined_bin = peak_idx + p

    freq = refined_bin * fs / N

    return freq


# ============================================================
# Parameters
# ============================================================

np.random.seed(42)

fs = 20000
duration = 5.0

t = np.arange(0, duration, 1/fs)

# ============================================================
# Frequency modulation
# ============================================================

f0 = 90.0

freq = np.ones_like(t) * 90.0
#freq = f0 + 2.0 * np.sin(2 * np.pi * 0.5 * t)

phase = 2 * np.pi * np.cumsum(freq) / fs

signal = np.sin(phase)

# ============================================================
# Add noise
# ============================================================

noisy = signal.copy()

# Gaussian noise
noisy += 0.3 * np.random.randn(len(noisy))

# Impulsive noise
num_spikes = 3200

idx = np.random.randint(0, len(noisy), num_spikes)

noisy[idx] += np.random.randn(num_spikes) * 8.0

# ============================================================
# Manual STFT
# ============================================================

window_size = 2048
hop_size = 128

window = np.hanning(window_size)
window_enh = np.hanning(int(window_size/2))

freqs = np.fft.rfftfreq(window_size, 1/fs)
freqs_enh = np.fft.rfftfreq(int(window_size/2), 1/fs)

num_frames = (len(noisy) - window_size) // hop_size

# STFT matrices
stft_raw = []
stft_enhanced = []

F0_ests_raw = []
F0_ests_enh = []

time_axis = []

for frame in range(num_frames):

    start = frame * hop_size
    end = start + window_size
    calc_data = noisy[start:end].copy()

    calc_data_posi = np.zeros_like(calc_data)
    calc_data_nega = np.zeros_like(calc_data)
    calc_data_posi[calc_data > 0] = calc_data[calc_data > 0]
    calc_data_nega[calc_data < 0] = -calc_data[calc_data < 0]
    segment_raw = calc_data * window
    periodicity_value = periodicity(calc_data_posi)+ periodicity(calc_data_nega)
    segment_enh = (periodicity_value - np.mean(periodicity_value)) * window_enh

    spec_raw = np.abs(np.fft.rfft(segment_raw))
    spec_enh = np.abs(np.fft.rfft(segment_enh))
    F0_raw = fft_peak_frequency(calc_data,fs=fs,fmin=70,fmax=110)
    F0_enh = fft_peak_frequency(periodicity_value,fs=fs,fmin=70,fmax=110)

    F0_ests_raw.append(F0_raw)
    F0_ests_enh.append(F0_enh)

    #plt.figure()
    #plt.plot(freqs[freqs < 200],spec_raw[freqs < 200])
    #plt.plot(freqs_enh[freqs_enh < 200],spec_enh[freqs_enh < 200])
    #plt.show()

    stft_raw.append(spec_raw)
    stft_enhanced.append(spec_enh)

    center_time = (start + window_size // 2) / fs
    time_axis.append(center_time)

stft_raw = np.array(stft_raw).T
stft_enhanced = np.array(stft_enhanced).T

time_axis = np.array(time_axis)

# dB scale
stft_raw_db = 20 * np.log10(stft_raw + 1e-12)
stft_enhanced_db = 20 * np.log10(stft_enhanced + 1e-12)

bedcmm_pitch,_ = bedcmmPitch.calc_Pitch(noisy,
                                        fs=fs,
                                        window_size=window_size,
                                        fmin=70,
                                        fmax=110,
                                        hop_size=hop_size,
                                        bedcmm_smooth=11)
                                        

# ============================================================
# Plot
# ============================================================

plt.figure(figsize=(10, 6))

# ------------------------------------------------------------
# Raw waveform
# ------------------------------------------------------------

plt.subplot(2, 3, 1)

plt.plot(t[:4000], noisy[:4000])

plt.title("Noisy signal")

plt.xlabel("Time [s]")

# ------------------------------------------------------------
# STFT raw
# ------------------------------------------------------------

plt.subplot(2, 3, 2)

plt.pcolormesh(
    time_axis,
    freqs,
    stft_raw_db,
    shading='gouraud'
)

plt.ylim(0, 200)

plt.title("Manual STFT - Raw")

plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")

plt.colorbar(label="dB")

# ------------------------------------------------------------
# STFT enhanced
# ------------------------------------------------------------

plt.subplot(2, 3, 5)

plt.pcolormesh(
    time_axis,
    freqs_enh,
    stft_enhanced_db,
    shading='gouraud'
)

plt.ylim(0, 200)

plt.title("Manual STFT - Periodicity extracted")

plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")

plt.colorbar(label="dB")

# ------------------------------------------------------------
# F0 est raw
# ------------------------------------------------------------

plt.subplot(2, 3, 3)

plt.plot(
    time_axis,
    F0_ests_raw,
    label=f"mean:{np.mean(F0_ests_raw):.3f},std:{np.std(F0_ests_raw):.3f}"
)

plt.ylim(85, 95)

plt.title("Manual F0 est - raw")

plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")

plt.legend()

# ------------------------------------------------------------
# STFT enhanced
# ------------------------------------------------------------

plt.subplot(2, 3, 6)

plt.plot(
    time_axis,
    F0_ests_enh,
    label=f"mean:{np.mean(F0_ests_enh):.3f},std:{np.std(F0_ests_enh):.3f}"
)

plt.ylim(85, 95)

plt.title("Manual F0 est - Periodicity extracted")

plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.legend()
# ------------------------------------------------------------
# F0 est bedcmmPitch
# ------------------------------------------------------------

plt.subplot(2, 3, 4)

plt.plot(
    time_axis,
    bedcmm_pitch[:-1],
    label=f"mean:{np.mean(bedcmm_pitch):.3f},std:{np.std(bedcmm_pitch):.3f}"
)

plt.ylim(85, 95)

plt.title("Manual F0 est - bedcmmPitch")

plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")

plt.legend()

plt.tight_layout()

plt.show()