心电信号stockwell变换代码

以下是心电信号进行Stockwell变换的Python代码示例：

import numpy as np
import pywt
import matplotlib.pyplot as plt

def st(x):
    # 计算Stockwell变换
    N = len(x)
    freqs = np.arange(1, N + 1) / float(N)
    tapers = np.fft.fft(np.exp(-1j * 2 * np.pi * np.arange(N) / N) / np.sqrt(N))
    tapers = np.outer(tapers, np.ones(N))
    X = np.fft.fft(x * tapers)
    S = np.zeros((N, N), dtype=np.complex128)
    for ii in range(N):
        S[:, ii] = pywt.scaleogram.cwt_frequencies(freqs, wavelet='cmor1.5-1', scales=[freqs[ii]])[0]
    S = np.abs(X)**2 / S / N
    return S

# 生成ECG信号
fs = 1000
x = np.zeros(int(5*fs))
x[int(0.5*fs):int(1.0*fs)] = 1.0
x[int(1.5*fs):int(2.0*fs)] = 0.5
x[int(2.5*fs):int(3.0*fs)] = 2.0
x[int(3.5*fs):int(4.0*fs)] = 1.5

# 计算Stockwell变换
S = st(x)

# 绘制结果
fig, (ax0, ax1) = plt.subplots(nrows=2, figsize=(6, 6))
im0 = ax0.plot(x)
ax0.set_title('ECG Signal')
ax0.set_xlabel('Time (ms)')
ax0.set_ylabel('Amplitude')
im1 = ax1.imshow(S, extent=[0, 5, 0, 500], aspect='auto', cmap='jet')
ax1.set_title('Stockwell Transform')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Frequency (Hz)')
fig.colorbar(im1, ax=ax1)
plt.show()

这里使用了PyWavelets库中的连续小波变换（CWT）计算不同尺度下的频率，然后根据Stockwell变换的定义计算每个时刻的时间-频率分布。最后使用matplotlib库绘制ECG信号和Stockwell变换结果。