对通道AI1-03的滤波数据进行希尔伯特黄变换,并将计算的总能量新建一列存放在"C:\Users\FILWYZ\Desktop\3.xlsx中。该表格中不需要保存“原始加速度第一个值,滤波后加速度第一个值”这两列数据,希尔伯特黄变换借鉴以下的代码
# ========== 希尔伯特黄变换核心代码 ==========
def find_extrema_matlab(x):
dx = np.diff(x)
d1 = dx[:-1]
d2 = dx[1:]
indmin = np.where((d1 * d2 < 0) & (d1 < 0))[0] + 1
indmax = np.where((d1 * d2 < 0) & (d1 > 0))[0] + 1
if np.any(dx == 0):
imax, imin = [], []
bad = (dx == 0)
dd = np.diff(np.concatenate([[0], bad.astype(int), [0]]))
debs = np.where(dd == 1)[0]
fins = np.where(dd == -1)[0]
for deb, fin in zip(debs, fins):
mid = int(np.round((deb + fin) / 2))
if dx[deb-1] > 0 and dx[fin-1] < 0:
imax.append(mid)
elif dx[deb-1] < 0 and dx[fin-1] > 0:
imin.append(mid)
if len(imax) > 0: indmax = np.sort(np.concatenate([indmax, imax]))
if len(imin) > 0: indmin = np.sort(np.concatenate([indmin, imin]))
return indmin, indmax
def boundary_conditions_matlab(indmin, indmax, t, x, nbsym=2):
lx = len(x)
if len(indmin) + len(indmax) < 3: raise Exception("not enough extrema")
if indmax[0] < indmin[0]:
if x[0] > x[indmin[0]]:
lmax = indmax[1: min(len(indmax), nbsym+1)][::-1]
lmin = indmin[:min(len(indmin), nbsym)][::-1]
lsym = indmax[0]
else:
lmax = indmax[:min(len(indmax), nbsym)][::-1]
lmin = np.concatenate((indmin[:min(len(indmin), nbsym-1)][::-1], [0]))
lsym = 0
else:
if x[0] < x[indmax[0]]:
lmax = indmax[:min(len(indmax), nbsym)][::-1]
lmin = indmin[1:min(len(indmin), nbsym+1)][::-1]
lsym = indmin[0]
else:
lmax = np.concatenate((indmax[:min(len(indmax), nbsym-1)][::-1], [0]))
lmin = indmin[:min(len(indmin), nbsym)][::-1]
lsym = 0
if indmax[-1] < indmin[-1]:
if x[-1] < x[indmax[-1]]:
rmax = indmax[max(len(indmax)-nbsym, 0):][::-1]
rmin = indmin[max(len(indmin)-nbsym, 0):-1][::-1]
rsym = indmin[-1]
else:
rmax = np.concatenate(([lx-1], indmax[max(len(indmax)-nbsym+1, 0):][::-1]))
rmin = indmin[max(len(indmin)-nbsym+1, 0):][::-1]
rsym = lx-1
else:
if x[-1] > x[indmin[-1]]:
rmax = indmax[max(len(indmax)-nbsym, 0):-1][::-1]
rmin = indmin[max(len(indmin)-nbsym+1, 0):][::-1]
rsym = indmax[-1]
else:
rmax = indmax[max(len(indmax)-nbsym+1, 0):][::-1]
rmin = np.concatenate(([lx-1], indmin[max(len(indmin)-nbsym+1, 0):][::-1]))
rsym = lx-1
tlmin = 2*t[lsym] - t[lmin]
tlmax = 2*t[lsym] - t[lmax]
trmin = 2*t[rsym] - t[rmin]
trmax = 2*t[rsym] - t[rmax]
zlmax = x[lmax]; zlmin = x[lmin]; zrmax = x[rmax]; zrmin = x[rmin]
tmin = np.concatenate((tlmin, t[indmin], trmin))
tmax = np.concatenate((tlmax, t[indmax], trmax))
zmin = np.concatenate((zlmin, x[indmin], zrmin))
zmax = np.concatenate((zlmax, x[indmax], zrmax))
return tmin, tmax, zmin, zmax
def interp_env_matlab(t_extrema, x_extrema, t_all, kind='spline'):
order = np.argsort(t_extrema)
t_extrema = t_extrema[order]; x_extrema = x_extrema[order]
_, unique_indices = np.unique(t_extrema, return_index=True)
t_extrema = t_extrema[unique_indices]; x_extrema = x_extrema[unique_indices]
if len(t_extrema) < 2: return np.full_like(t_all, x_extrema[0] if len(x_extrema) > 0 else 0)
if kind == 'spline':
if len(t_extrema) >= 4:
cs = CubicSpline(t_extrema, x_extrema, extrapolate=True)
return cs(t_all)
else:
f = interp1d(t_extrema, x_extrema, kind='linear', fill_value='extrapolate')
return f(t_all)
elif kind == 'pchip':
if len(t_extrema) >= 2:
pchip = PchipInterpolator(t_extrema, x_extrema, extrapolate=True)
return pchip(t_all)
else:
return np.full_like(t_all, np.mean(x_extrema))
elif kind == 'linear':
f = interp1d(t_extrema, x_extrema, kind='linear', fill_value='extrapolate')
return f(t_all)
else:
raise ValueError("Unknown interp kind")
def emd_m_full(x, t=None, stop=[0.05,0.5,0.05], maxiter=2000, maxmodes=None, interp='spline', nbsym=2):
if t is None: t = np.arange(len(x))
x = np.asarray(x, dtype=np.float64); N = len(x); imfs = []; r = x.copy()
sd, sd2, tol = stop; k = 0; IMF_limit = maxmodes if maxmodes is not None else N
while True:
indmin, indmax = find_extrema_matlab(r)
if len(indmin) + len(indmax) < 3 or (IMF_limit > 0 and k >= IMF_limit): break
m = r.copy()
for _ in range(maxiter):
indmin, indmax = find_extrema_matlab(m)
if len(indmin) + len(indmax) < 3: break
tmin, tmax, zmin, zmax = boundary_conditions_matlab(indmin, indmax, t, m, nbsym=nbsym)
envmin = interp_env_matlab(tmin, zmin, t, kind=interp)
envmax = interp_env_matlab(tmax, zmax, t, kind=interp)
mean_env = (envmin + envmax) / 2
amp = np.abs(envmax - envmin) / 2
mean_amp = np.abs(mean_env)
sx = mean_amp / (amp + 1e-12)
stop1 = (np.mean(sx > sd) > tol)
stop2 = np.any(sx > sd2)
stop3 = np.abs(len(indmin) - len(indmax)) > 1
if not ((stop1 or stop2) and not stop3): break
m = m - mean_env
imfs.append(m)
r = r - m
k += 1
if np.any(np.abs(r) > 1e-10): imfs.append(r)
return np.array(imfs)
def hhspectrum(imfs, t, l=1, dt=1/1000):
Nmodes, N = imfs.shape
tt = t[l:N-l]
A = np.zeros((Nmodes, len(tt)))
F = np.zeros((Nmodes, len(tt)))
for i in range(Nmodes):
analytic_signal = hilbert(imfs[i, :])
amplitude_envelope = np.abs(analytic_signal)[l:N-l]
phase = np.unwrap(np.angle(analytic_signal))
freq = (phase[l+1:N] - phase[l-1:N-2]) / (4.0 * np.pi * dt)
A[i, :] = amplitude_envelope
F[i, :] = freq
return A, F, tt