-
Notifications
You must be signed in to change notification settings - Fork 0
/
analysis.py
138 lines (99 loc) · 2.76 KB
/
analysis.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
"""
File: analysis.py
Author: Chuncheng Zhang
Date: 2024-03-04
Copyright & Email: chuncheng.zhang@ia.ac.cn
Purpose:
Analysis the signal loaded from load_data.py
Functions:
1. Requirements and constants
2. Function and class
3. Play ground
4. Pending
5. Pending
"""
# %% ---- 2024-03-04 ------------------------
# Requirements and constants
import pandas as pd
import plotly.express as px
import matplotlib.pyplot as plt
import numpy as np
from scipy import fft
from scipy import signal
from pathlib import Path
from load_data import load_data, fs, ch_names
# %% ---- 2024-03-04 ------------------------
# Function and class
window_names = [
e.strip()
for e in """
boxcar
triang
blackman
hamming
hann
bartlett
flattop
parzen
bohman
blackmanharris
nuttall
barthann
cosine
exponential
tukey
taylor
lanczos
""".split(
"\n"
)
if e.strip()
]
def compute_psd(window: str = "hann", **kwargs):
data = load_data(Path("data/10hz/S0.mat"))
print(kwargs)
f, Pxx_den = signal.welch(data, fs, window, **kwargs)
df = mk_df(Pxx_den, ch_names, f)
return f, Pxx_den, df
def mk_df(Pxx_den, ch_names, f):
array = []
for i in range(Pxx_den.shape[0]):
for j, ch in enumerate(ch_names):
for k, freq in enumerate(f):
array.append(
dict(epoch_idx=i, ch_name=ch, freq=freq, value=Pxx_den[i][j][k])
)
df = pd.DataFrame(array)
return df
# %% ---- 2024-03-04 ------------------------
# Play ground
if __name__ == "__main__":
data = load_data(Path("data/4hz/S0.mat"))
# --------------------
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.plot(data[0][0])
# --------------------
f_1, Pxx_den_1 = signal.welch(data, fs, window="hamming")
f_2, Pxx_den_2 = signal.welch(data, fs, window="hann")
diff = np.abs(Pxx_den_1 - Pxx_den_2)
f = f_1
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
for j in range(diff.shape[0]):
ax.semilogy(f, diff[j][0], label=f"Epochs {j}")
ax.set_xlabel("frequency [Hz]")
ax.set_ylabel("PSD [V**2/Hz]")
ax.set_title(f"Difference between hamming and hann window")
# --------------------
f, Pxx_den = signal.welch(data, fs, window="hamming")
print(Pxx_den.shape)
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
for j in range(Pxx_den.shape[0]):
ax.semilogy(f, Pxx_den[j][0], label=f"Epochs {j}")
ax.set_xlabel("frequency [Hz]")
ax.set_ylabel("PSD [V**2/Hz]")
plt.show()
# %%
# %% ---- 2024-03-04 ------------------------
# Pending
# %% ---- 2024-03-04 ------------------------
# Pending