Python 过滤 4D 结构中的信号

Question

我有一个 4D 结构，其中包含形状为 (3432,1,30,512) 的 EEG 数据

这代表 3432 次数据试验，其中每个试验包含 30 个电极的 512 个时间样本

对于每次试验，我想用某些过滤参数过滤每个电极，从而产生一个新的 2D (30,512) 过滤数组。然后，我想将它们沿 4D 结构的轴 1 堆叠。

目前，我正在迭代每个试验，然后是每个电极，过滤信号并将所有内容重新插入 4D：

from scipy import signal    

NUM_CHANS = 30
NUM_TIMESTAMPS = 512
FREQ_BANDS = ((0.1, 3), (4, 7), (8, 13), (16, 31))

all_data = np.reshape(all_data, (-1, 1, NUM_CHANS, NUM_TIMESTAMPS))  # (3432,1,30,512)

for i in range(all_data.shape[0]):
    num_band = 1
    for band in FREQ_BANDS:
        lower = float(band[0])/(SAMPLE_RATE/2)
        upper = float(band[1])/(SAMPLE_RATE/2)

        # Design new filter for the current frequency band
        b, a = signal.butter(2, [lower, upper], 'bandpass')

        temp_trial = np.zeros((NUM_CHANS, NUM_TIMESTAMPS))

        for ch in range(NUM_CHANS):
            # Filter the current electrode
            output_signal = signal.filtfilt(b, a, all_data[i,0,ch,:])
            temp_trial[ch,:] = output_signal

        # Insert temp_trial (2D) into all_data (4D) along axis 1

        num_band += 1

迭代试验和电极非常慢（完成整个循环大约需要 2 小时）。是否有更有效的方法将此过滤器应用于所有电极/试验？

我一直在尝试找到一种在 2D 上应用过滤器的方法，因此我不需要迭代电极，但找不到任何东西。

编辑：

这是使用filtfilt 轴参数的正确方法吗？

all_data = np.reshape(all_data, (-1, 1, NUM_CHANS, NUM_TIMESTAMPS))

filtered_data = [all_data]

for band in FREQ_BANDS:
    lower = float(band[0])/(SAMPLE_RATE/2)
    upper = float(band[1])/(SAMPLE_RATE/2)
    b, a = signal.butter(2, [lower, upper], 'bandpass')

    output_signal = signal.filtfilt(b, a, all_data, axis=3)

    filtered_data.append(output_signal)

all_data = np.concatenate(filtered_data, axis=1)

Answer 1

这是一项可能会提高性能的更改。 filtfilt 有一个轴参数，它允许您沿 n 维数组的轴应用相同的一维滤波器。您可以替换此代码

    temp_trial = np.zeros((NUM_CHANS, NUM_TIMESTAMPS))
    for ch in range(NUM_CHANS):
        # Filter the current electrode
        output_signal = signal.filtfilt(b, a, all_data[i,0,ch,:])
        temp_trial[ch,:] = output_signal

与

    temp_trial = signal.filtfilt(b, a, all_data[i,0,:,:], axis=1)

由于默认的axis 是axis=-1，如果您愿意，可以省略轴参数：

    temp_trial = signal.filtfilt(b, a, all_data[i,0,:,:])

---

通过重新排列外部循环，您可以做得更好。使“带”环成为外环。然后，对于每个波段，您可以将 filtfilt once 应用于 3-d 数组 all_data[:, 0, :, :]。

类似这样的：

shp = all_data.shape

# This assumes `all_data` has shape (3432,1,30,512), but I don't
# think you need that trivial extra dimension in there if you
# preallocate `filtered_data` like I am doing here.
filtered_data = np.empty(shp[0:1] + (len(FREQ_BANDS),) + shp[2:])

for k, band in enumerate(FREQ_BANDS):
    lower = float(band[0])/(SAMPLE_RATE/2)
    upper = float(band[1])/(SAMPLE_RATE/2)

    # Design new filter for the current frequency band
    b, a = signal.butter(2, [lower, upper], 'bandpass')

    filtered_data[:, k, :, :] = filtfilt(b, a, all_data[:,0,:,:])

这不包括filtered_data 中的原始all_data。如果需要，请在创建filtered_data 的np.empty() 调用中使用len(FREQ_BANDS)+1，设置filtered_data[:,0,:,:] = all_data，并在保存filtfilt() 调用结果的赋值语句中使用k+1 而不是k .