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()
この記事が気に入ったらサポートをしてみませんか?