hilbert_example.py

#! /usr/bin/env python3

""" A Python script for demostration of Hilbert transform."""

# Copyright (C) 2020 by Akira TAMAMORI

# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

from scipy.io import wavfile
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal

WAVE_FILE = "./wavefile.wav"


def main():
    """Main function."""

    sample_rate, wave_data = wavfile.read(WAVE_FILE)

    # ヒルベルト変換 (FFT -> 虚部0 & 実部2倍 -> 逆FFT)
    envelop = np.abs(signal.hilbert(wave_data))  # 包絡
    angle = np.unwrap(np.angle(signal.hilbert(wave_data)))  # 瞬時位相
    reconst = envelop * np.cos(angle)  # 再構成

    # plot original and reconstructed waveform
    time = np.arange(0, len(wave_data)) / sample_rate
    fig = plt.figure(figsize=(10, 6))
    axes1 = fig.add_subplot(3, 1, 1)
    axes1.set_xlabel("Time [s]")
    axes1.set_ylabel("Amplitude")
    axes1.set_title("Original waveform")
    axes1.plot(time, wave_data)

    axes2 = fig.add_subplot(3, 1, 2)
    axes2.set_xlabel("Time [s]")
    axes2.set_ylabel("Amplitude")
    axes2.set_title("Reconstructed waveform")
    axes2.plot(time, reconst)

    axes3 = fig.add_subplot(3, 1, 3)
    axes3.set_xlabel("Time [s]")
    axes3.set_ylabel("Amplitude")
    axes3.set_title("Error")
    axes3.plot(time, reconst - wave_data)
    fig.tight_layout()
    plt.show()


if __name__ in '__main__':
    main()