Python Filtern von Signalen in einer 4D-Struktur

Question

Ich habe eine 4D-Struktur mit Form EEG-Daten enthält (3432,1,30,512)

Dies entspricht 3432 Studien von Daten, wobei jeder Versuch 512 Zeitproben für 30 Elektroden

Für jeden Versuch enthält, Ich möchte jede Elektrode mit bestimmten Filterparametern filtern, was zu einem neuen gefilterten Array führt. Dann möchte ich diese entlang der Achse 1 der 4D-Struktur stapeln.

Derzeit bin ich über jeden Versuch laufen, dann jede Elektrode, Filtern des Signals und alles wieder in die 4D Einfügen:

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 

Iterieren über Studien und Elektroden ist extrem langsam (dauert ca. 2 Stunden vervollständige die ganze Schleife). Gibt es einen effizienteren Weg diesen Filter über alle Elektroden/Versuche anzuwenden?

Ich habe versucht, einen Weg zu finden, den Filter über 2D anzuwenden, so dass ich nicht über Elektroden iterieren muss, aber nichts gefunden habe.

EDIT:

Wäre dies der richtige Weg sein filtfilt Achse Argument zu benutzen?

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

Hier ist eine Änderung, die die Leistung verbessern könnte. filtfilt hat ein Achsenargument, mit dem Sie den gleichen 1-d-Filter entlang einer Achse eines n-dimensionalen Arrays anwenden können. Sie können diesen Code

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 

mit

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

ersetzen Und da der Standard axisaxis=-1 ist, können Sie die Achse Argument auslassen, wenn Sie bevorzugen:

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

---

Sie kann noch besser, indem Sie die äußeren Schleifen neu anordnen. Machen Sie die "Band" -Schleife zur äußeren Schleife. Dann können Sie für jedes Band filtfilteinmal auf das 3-d-Array anwenden.

Etwas wie folgt aus:

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,:,:]) 

Dies gilt nicht für die ursprüngliche all_data in filtered_data. Wenn Sie das möchten, verwenden Sie len(FREQ_BANDS)+1 in dem Aufruf an np.empty(), das filtered_data erstellt, filtered_data[:,0,:,:] = all_data setzt und k+1 anstelle von k in Zuordnungsanweisung verwendet, die das Ergebnis des filtfilt()-Aufrufs speichert.