4D结构中的Python滤波信号

Question

我有一个包含(3432,1,30,512)形状的脑电数据的4D结构

这代表了3432次数据试验，每次试验包含512次30个电极的样本。

在每一次试验中，我都希望用一定的滤波参数过滤每个电极，这就产生了一个新的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有一个axis参数，它允许您沿着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，所以如果您愿意的话，可以省略axis参数：

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

通过重新排列外循环，你可以做得更好。将“带”环设为外环。然后，对于每个波段，您可以将filtfilt一次应用于三维数组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中的原始filtered_data。如果需要，在创建len(FREQ_BANDS)+1、设置filtered_data[:,0,:,:] = all_data和在赋值语句中使用k+1而不是k来保存filtfilt()调用结果的np.empty()调用中使用k+1。