python勉強

・seabornで綺麗にグラフを書く方法　sns.set()
・FFT
・Low Pass Filter

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal, interpolate
import seaborn as sns
import matplotlib as mpl

sns.set()
# mpl.style.use("ggplot")

# 移動平均
def LPF_MAM(x, k):
    x_mean = np.zeros(x.shape)
    N = x.shape[0]
    for i in range(N):
        if i - k // 2 < 0:
            x_mean[i] = x[: i - k // 2 + k].mean()
        elif i - k // 2 + k >= N:
            x_mean[i] = x[i - k // 2 :].mean()
        else:
            x_mean[i] = x[i - k // 2 : i - k // 2 + k].mean()
    return x_mean


def LPF_CF(x, times, fmax):
    freq_X = np.fft.fftfreq(times.shape[0], times[1] - times[0])
    X_F = np.fft.fft(x)
    X_F[freq_X > fmax] = 0
    X_F[freq_X < -fmax] = 0
    #     虚数は削除
    x_CF = np.fft.ifft(X_F).real
    return x_CF


def resample(x, y, num):
    f = interpolate.interp1d(x, y, kind="cubic")
    x_resample = np.linspace(x[0], x[-1], num)
    y_resample = f(x_resample)
    return x_resample, y_resample


# FFT
def fft_(dx, y):
    freq = np.linspace(0, int(1.0 / dx), int(y.shape[0]))  # 周波数軸
    F = np.fft.fft(y)
    # 正規化 + 交流成分2倍
    F = F / (y.shape[0] / 2)
    F[0] = F[0] / 2
    return freq, F


def LPF_GC(x, times, sigma):
    sigma_k = sigma / (times[1] - times[0])
    kernel = np.zeros(int(round(3 * sigma_k)) * 2 + 1)
    for i in range(kernel.shape[0]):
        kernel[i] = (
            1.0
            / np.sqrt(2 * np.pi)
            / sigma_k
            * np.exp((i - round(3 * sigma_k)) ** 2 / (-2 * sigma_k ** 2))
        )

    kernel = kernel / kernel.sum()
    x_long = np.zeros(x.shape[0] + kernel.shape[0])
    x_long[kernel.shape[0] // 2 : -kernel.shape[0] // 2] = x
    x_long[: kernel.shape[0] // 2] = x[0]
    x_long[-kernel.shape[0] // 2 :] = x[-1]

    x_GC = np.convolve(x_long, kernel, "same")

    return x_GC[kernel.shape[0] // 2 : -kernel.shape[0] // 2]


y = np.array(
    [
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        5,
        3,
        20,
        15,
        30,
        50,
        60,
        65,
        59,
        62,
        40,
        30,
        33,
        20,
        10,
        5,
        6,
        4,
        2,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
    ]
)


dx = 0.5
N = y.shape[0] / 2
print(y.shape[0])
x = np.arange(0, N, dx)

fig, ax = plt.subplots()
ax.plot(x, y, alpha=0.7, linewidth=2)
ax.plot(x, LPF_MAM(y, 3), alpha=0.7, linewidth=2)
ax.plot(x, LPF_CF(y, x, 0.25), alpha=0.7, linewidth=2)
# ax.scatter(x, y, alpha=0.7)
# ax.scatter(x, LPF_MAM(y, 3), alpha=0.7)
# ax.scatter(x, LPF_CF(y, x, 0.25), alpha=0.7)
plt.show()

fig, ax = plt.subplots()
ax.plot(fft_(dx, y)[0], np.abs(fft_(dx, y)[1]), alpha=0.7)
ax.scatter(fft_(dx, y)[0], np.abs(fft_(dx, y)[1]), alpha=0.7)
plt.show()

# リサンプリング
num = 100
x_resampled = resample(x, y, num)[0]
y_resampled = resample(x, y, num)[1]
fig, ax = plt.subplots()
ax.plot(x, y, alpha=0.7, linewidth=2)
ax.plot(x_resampled, y_resampled, alpha=0.7, linewidth=2)
plt.show()

freq, F = fft_(x_resampled[1] - x_resampled[0], y_resampled)
fig, ax = plt.subplots()
ax.plot(freq, np.abs(F), alpha=0.7)
ax.scatter(freq, np.abs(F), alpha=0.7)
plt.show()

fig, ax = plt.subplots()
ax.plot(x_resampled, y_resampled, alpha=0.7, linewidth=2)
ax.plot(x_resampled, LPF_MAM(y_resampled, 10), alpha=0.7, linewidth=2)
ax.plot(x_resampled, LPF_GC(y_resampled, x_resampled, 2), alpha=0.7, linewidth=2)
# ax.scatter(x, y, alpha=0.7)
# ax.scatter(x, LPF_MAM(y, 3), alpha=0.7)
# ax.scatter(x, LPF_CF(y, x, 0.25), alpha=0.7)
plt.show()

exit()