add summary script

.gitignore

@@ -205,3 +205,6 @@ cython_debug/
 marimo/_static/
 marimo/_lsp/
 __marimo__/
+
+# wykresy danych
+data/plots/**

data/loader.py

@@ -0,0 +1,65 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import mne
+
+
+def load_raw_data(
+    filename: str, l_freq: float = 6.0, h_freq: float = 30.0
+) -> mne.io.Raw:
+    try:
+        raw = mne.io.read_raw_fif(filename, preload=True)
+        raw = raw.load_data()
+        raw = raw.pick_types(eeg=True, stim=False, eog=False, exclude="bads")
+        raw.apply_function(lambda x: x * 10**-6)
+        raw.filter(l_freq=l_freq, h_freq=h_freq)
+        raw.notch_filter(freqs=50)
+        return raw
+    except FileNotFoundError:
+        print(f"Nie można znaleźć pliku '{filename}'.")
+        return None
+
+
+def load_data(
+    filename: str,
+    target_marker_freqs: list[float],
+    l_freq: float = 6.0,
+    h_freq: float = 30.0,
+    window_time_frame: int = 5,
+) -> dict[np.ndarray] | None:
+    raw = load_raw_data(filename, l_freq=l_freq, h_freq=h_freq)
+    results = {}
+    if raw is not None:
+        events, event_dict = mne.events_from_annotations(raw)
+        print("Dostępne markery w nagraniu:", event_dict)
+
+        for target_marker_freq in target_marker_freqs:
+            target_event_id = None
+            for marker_name, marker_id in event_dict.items():
+                if f"TARGET_{target_marker_freq}" in marker_name:
+                    target_event_id = marker_id
+                    print(f"Wybrano marker: {marker_name} (ID: {marker_id})")
+                    break
+
+            if target_event_id is not None:
+                epochs = mne.Epochs(
+                    raw,
+                    events,
+                    event_id=target_event_id,
+                    tmin=0.0,
+                    tmax=window_time_frame,
+                    baseline=None,
+                    preload=True,
+                )
+
+                target_epoch_data = epochs.get_data()[0]
+                print(f"Wyizolowano epoch o wymiarach: {target_epoch_data.shape}")
+                results[target_marker_freq] = target_epoch_data
+            else:
+                print(
+                    f"Nie odnaleziono markera dla bodźca: {target_marker_freq} Hz w pliku!"
+                )
+                results[target_marker_freq] = raw.get_data()
+    else:
+        return None
+
+    return results

data/plots.py

@@ -0,0 +1,81 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fft import fft, fftfreq
+
+plt.rcParams["figure.figsize"] = (12, 6)
+plt.style.use("seaborn-v0_8-darkgrid")
+
+
+def plot_time_series(
+    time, raw_signal, car_signal, filtered_signal, ch_idx, sg_window=21, sg_polyorder=3
+):
+    plt.figure(figsize=(12, 6))
+
+    plt.plot(time, raw_signal, label="Sygnał Surowy (Raw)", alpha=0.5, color="gray")
+    plt.plot(time, car_signal, label="Sygnał po CAR", alpha=0.5, color="blue")
+    plt.plot(
+        time,
+        filtered_signal,
+        label=f"SG Filter (win={sg_window}, order={sg_polyorder})",
+        linewidth=1.5,
+        color="red",
+    )
+
+    plt.title(f"Cała wyodrębniona epoka stymulacji (Kanał {ch_idx})")
+    plt.xlabel("Czas względny stymulacji [s]")
+    plt.ylabel("Amplituda [V]")
+    plt.legend()
+    plt.tight_layout()
+
+
+def plot_fft(signal, title, color, sampling_rate=250):
+    # Obliczanie widma amplitudowego
+    N = len(signal)
+    signal_fft = fft(signal)
+    freqs = fftfreq(N, 1 / sampling_rate)
+
+    # Bierzemy pod uwagę tylko częstotliwości dodatnie
+    pos_mask = freqs > 0
+    freqs = freqs[pos_mask]
+    amplitudes = np.abs(signal_fft[pos_mask]) / N
+
+    # Zawężanie widma do SSVEP
+    mask_ssvep = freqs <= 40
+
+    plt.plot(
+        freqs[mask_ssvep], amplitudes[mask_ssvep], label=title, color=color, alpha=0.8
+    )
+    plt.xlabel("Częstotliwość (Hz)")
+    plt.ylabel("Amplituda")
+    plt.legend()
+    plt.tight_layout()
+
+
+def plot_psd(f_raw, psd_raw, f_filt, psd_filt):
+    # Widmo SSVEP max 40 Hz
+    plot_range_raw = f_raw <= 40
+    plot_range_filt = f_filt <= 40
+
+    plt.figure(figsize=(12, 6))
+    plt.semilogy(
+        f_raw[plot_range_raw],
+        psd_raw[plot_range_raw],
+        label="Surowy Sygnał",
+        color="blue",
+        alpha=0.5,
+    )
+    plt.semilogy(
+        f_filt[plot_range_filt],
+        psd_filt[plot_range_filt],
+        label=f"Wzbogacony Sygnał (SG Filter)",
+        color="red",
+        linewidth=1.5,
+    )
+
+    plt.title(
+        "Power Spectral Density (Spodziewane potężne 'piki' przy targetach SSVEP)"
+    )
+    plt.xlabel("Częstotliwość [Hz]")
+    plt.ylabel("PSD [V^2/Hz]")
+    plt.legend()
+    plt.tight_layout()

data/summary.py

@@ -0,0 +1,152 @@
+import sys
+import os
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+project_root = os.path.abspath(os.path.join(current_dir, ".."))
+if project_root not in sys.path:
+    sys.path.insert(0, project_root)
+from experiment.src.pipeline.predictor import CCAPredictor
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import savgol_filter, welch
+from sklearn.cross_decomposition import CCA
+from loader import load_data
+from plots import plot_time_series, plot_fft, plot_psd
+
+from experiment.src.config.bci_config import (
+    SAMPLING_RATE,
+    TARGET_CHANNELS,
+    SG_WINDOW,
+    SG_POLYORDER,
+    FREQS,
+    WINDOW_TIME_FRAME,
+)
+from experiment.src.pipeline.predictor import CCAPredictor
+
+def main():
+    filenames = [
+        "kacper/kacper1.fif",
+        "kacper/kacper2.fif",
+        "kacper/kacper3.fif",
+    ]
+
+    for filename in filenames:
+        print(f"\n{'='*50}\nPrzetwarzanie pliku: {filename}\n{'='*50}")
+        data_dict = load_data(filename, target_marker_freqs=FREQS, l_freq=5.0, h_freq=45.0)
+
+        if data_dict is not None:
+            # Iteracja po słowniku z wyizolowanymi epokami dla każdej częstotliwości
+            for target_freq, epoch_data in data_dict.items():
+                print(
+                    f"\n--- Analiza dla częstotliwości docelowej: {target_freq} Hz ---"
+                )
+
+                # Zabezpieczenie na wypadek, gdyby epoki nie pobrano prawidłowo
+                if epoch_data.ndim < 2:
+                    print(f"Pominięto {target_freq} Hz (brak danych).")
+                    continue
+
+                num_samples = epoch_data.shape[1]
+                time = np.linspace(0, WINDOW_TIME_FRAME, num_samples, endpoint=False)
+
+                # Operacja CAR (Common Average Reference) dla całej wyizolowanej epoki
+                mean_signal = np.mean(epoch_data[TARGET_CHANNELS, :], axis=0)
+                car_all_channels = epoch_data - mean_signal
+
+                # Analiza per kanał (Time Series, FFT, PSD)
+                for ch_idx in TARGET_CHANNELS:
+                    print(f"  -> Rysowanie wykresów dla kanału nr {ch_idx}")
+
+                    raw_signal = epoch_data[ch_idx, :]
+                    car_signal = car_all_channels[ch_idx, :]
+                    filtered_signal = savgol_filter(
+                        car_signal, window_length=SG_WINDOW, polyorder=SG_POLYORDER
+                    )
+
+                    # 1. Wykres w dziedzinie czasu
+                    plot_time_series(
+                        time, raw_signal, car_signal, filtered_signal, ch_idx
+                    )
+                    filepath = os.path.join(
+                        "plots",
+                        "time_series",
+                        filename,
+                        str(ch_idx),
+                        f"{target_freq}.png",
+                    )
+                    os.makedirs(os.path.dirname(filepath), exist_ok=True)
+                    plt.savefig(filepath)
+                    plt.close()
+
+                    # 2. Wykres FFT (ponieważ plot_fft rysuje jedną linię, grupujemy je w jednej figurze)
+                    plt.figure(figsize=(12, 6))
+                    plot_fft(raw_signal, "Przed filtrem", "gray")
+                    plot_fft(car_signal, "Po CAR", "blue")
+                    plot_fft(filtered_signal, "Po filtrze Savitzky-Golay", "red")
+                    plt.title(
+                        f"Analiza częstotliwości FFT - Kanał {ch_idx} (Bodziec: {target_freq} Hz)"
+                    )
+                    filepath = os.path.join(
+                        "plots", "fft", filename, str(ch_idx), f"{target_freq}.png"
+                    )
+                    os.makedirs(os.path.dirname(filepath), exist_ok=True)
+                    plt.savefig(filepath)
+                    plt.close()
+
+                    # 3. Wykres PSD (Welch)
+                    nperseg = min(2 * SAMPLING_RATE, len(raw_signal))
+                    f_raw, psd_raw = welch(
+                        raw_signal, fs=SAMPLING_RATE, nperseg=nperseg
+                    )
+                    f_filt, psd_filt = welch(
+                        filtered_signal, fs=SAMPLING_RATE, nperseg=nperseg
+                    )
+                    plot_psd(f_raw, psd_raw, f_filt, psd_filt)
+                    filepath = os.path.join(
+                        "plots", "psd", filename, str(ch_idx), f"{target_freq}.png"
+                    )
+                    os.makedirs(os.path.dirname(filepath), exist_ok=True)
+                    plt.savefig(filepath)
+                    plt.close()
+
+                # 4. CCA dla całej epoki (wykorzystuje zbiór kanałów docelowych)
+                print("  -> Generowanie klasyfikacji CCA")
+                epoch_X = car_all_channels[TARGET_CHANNELS, :].T
+
+                predictor = CCAPredictor(
+                    expected_frequencies=FREQS,
+                    sampling_rate=SAMPLING_RATE,
+                    window_length=num_samples,
+                    threshold=0.5,
+                )
+                predicted_freq = predictor.predict(epoch_X)
+                print(f"     Predykcja CCA: {predicted_freq} Hz (dla {target_freq} Hz)")
+                print(
+                    f"     Współczynniki korelacji CCA: {predictor.last_correlations}"
+                )
+                freqs_list = list(predictor.last_correlations.keys())
+                corr_values = list(predictor.last_correlations.values())
+
+                plt.figure(figsize=(10, 5))
+                bars = plt.bar(
+                    [str(f) for f in freqs_list], corr_values, color="teal", alpha=0.7
+                )
+
+                best_idx = np.argmax(corr_values)
+                bars[best_idx].set_color("crimson")
+
+                plt.title(
+                    f"Klasyfikacja CCA dla epoki {target_freq} Hz (na podstawie {len(TARGET_CHANNELS)} kanałów)"
+                )
+                plt.xlabel("Referencyjne częstotliwości z bazy SSVEP [Hz]")
+                plt.ylabel("Współczynnik korelacji CCA")
+                plt.tight_layout()
+                filepath = os.path.join("plots", "cca", filename, f"{target_freq}.png")
+                os.makedirs(os.path.dirname(filepath), exist_ok=True)
+                plt.savefig(filepath)
+                plt.close()
+
+
+if __name__ == "__main__":
+    main()

experiment/src/pipeline/predictor.py

@@ -25,6 +25,8 @@ def __init__(
         self.reference_signals_dict = self._generate_reference_signals_dictionary(
             window_length
         )
+        self.last_correlations: Dict[float, float] = {}
+        self.max_correlation: Optional[float] = None
 
     def _generate_reference_signals_for_frequency(
         self, length: int, freq: float, num_harmonics=2
@@ -60,11 +62,14 @@ def predict(self, X_preprocessed: np.ndarray) -> Optional[float]:
             Xc, Yc = self.cca.transform(X_preprocessed, Y)
 
             corr = np.corrcoef(Xc[:, 0], Yc[:, 0])[0, 1]
+            self.last_correlations[expected_freq] = corr
 
             if corr > max_corr:
                 max_corr = corr
                 best_freq = expected_freq
 
+        self.max_correlation = max_corr
+
         if max_corr > self.threshold:
             return best_freq