Pythonで時刻歴データの包絡線を書く

時刻歴データを扱っている時，ピーク値を俯瞰したようなグラフをプロットしたい．

envelope_both.py

# -*- coding: utf-8 -*-
import numpy as np
from scipy import fftpack
from matplotlib import pyplot as plt
 
# FFTをする関数
def fft_ave(data, samplerate, Fs):
    fft = fftpack.fft(data)                      # FFT（実部と虚部）
    fft_amp = np.abs(fft / (Fs / 2))             # 振幅成分を計算
    fft_axis = np.linspace(0, samplerate, Fs)    # 周波数軸を作成
    return fft, fft_amp, fft_axis
 
# サンプルの時間領域信号を生成
dt = 0.001
t = np.arange(0, 1, dt)
wave = np.cos(2 * np.pi * 50 * t)
wave *= (1 + 0.5 * np.sin(2 * np.pi * 2 * t))
wave *= (1 + 0.5 * np.sin(2 * np.pi * 10 * t))
wave += 1
wave_ave = np.average(wave)
plt.figure()
plt.plot(t,wave)
plt.show()
wave = wave - wave_ave

# FFTする
fft, fft_amp, fft_axis = fft_ave(wave, 1 / dt, len(wave))
 
# 負周波数域（ナイキスト周波数以降）を0にする（実部と虚部の両方）
zeros = np.zeros(int(len(fft) / 2))
fft[int(len(fft) / 2):len(fft)] = zeros
 
# 正周波数域（ナイキスト周波数まで）を2倍する
fft *= 2
 
# 操作した実部と虚部を持つ周波数波形をIFFTし、絶対値をとる（包絡線を得る）
ifft_time = np.abs(fftpack.ifft(fft))
 
# 下側の包絡線を計算
lower_envelope = -ifft_time
 
# ここからグラフ描画
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'Times New Roman'
 
# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
 
# グラフの上下左右に目盛線を付ける。
plt.figure()
plt.plot(t,wave+wave_ave, label='original', lw=1)
plt.plot(t, ifft_time+wave_ave, label='upper envelope', lw=1)
plt.plot(t, lower_envelope+wave_ave, label='lower envelope', lw=1)

plt.legend()
 
# グラフを表示する。
plt.show()
plt.close()