spc_analyse.py

import os
import numpy as np
import matplotlib.pyplot as plt
import edfio
from matplotlib.lines import Line2D

path1 = '/home/data/data/seizeIT2/ds005873-1.0.1/sub-002/ses-01/eeg/sub-002_ses-01_task-szMonitoring_run-01_eeg.edf'
eeg_data = edfio.read_edf(path1)
signals = np.array([sig.data for sig in eeg_data.signals], dtype=np.float32)

window_size = 512
stride = 256
k = 5
anomalies_per_channel = {}

for ch in range(signals.shape[0]):
    anomalies = []
    signal = signals[ch]

    for start in range(0, len(signal) - window_size, stride):
        window = signal[start:start + window_size]
        mu = np.mean(window)
        sigma = np.std(window)
        upper = mu + k * sigma
        lower = mu - k * sigma

        if np.any((window > upper) | (window < lower)):
            anomalies.append((start, mu, sigma))

    anomalies_per_channel[ch] = anomalies

for ch, anomalies in anomalies_per_channel.items():
    print(f'{ch}: {len(anomalies)} anomalies detected')
    for anomaly in anomalies:
        print(f'  Start: {anomaly[0]}, Mean: {anomaly[1]:.2f}, Std: {anomaly[2]:.2f}')


def plot_anomalies_time_series(time_series, anomalies, window_size, stride, save_path="graphs/anomalies_time_series", duration_plot_seconds=60, pre_seconds=5, center_seconds=None, start_sample=None):
    """
    Plot time series with detected anomalies, showing start and end of each anomaly window.

    Parameters:
    - time_series (numpy.ndarray): Shape (channels, samples)
    - anomalies (dict): Dict of channel -> list of (start_index, mean, std)
    - window_size (int): Size of the anomaly detection window (in samples)
    - stride (int): Stride used in anomaly detection (in samples)
    - save_path (str): Path to save the plot
    - duration_plot_seconds (int): Duration to plot in seconds
    - pre_seconds (int): Seconds before the first anomaly to start plotting (ignored if center_seconds is set)
    - center_seconds (float or None): If set, center the plot window around this second (ignored if start_sample is set)
    - start_sample (int or None): If set, use this sample index as the plot start

    Returns:
    int: The sample index of the first detected anomaly (or None if no anomaly found)
    """
    sampling_rate = 256
    first_anomaly = None
    for ch_anoms in anomalies.values():
        if ch_anoms:
            idx = min(a[0] for a in ch_anoms)
            if first_anomaly is None or idx < first_anomaly:
                first_anomaly = idx

    if start_sample is not None:
        start_cond = max(0, start_sample)
    elif center_seconds is not None:
        center_sample = int(center_seconds * sampling_rate)
        start_cond = max(0, center_sample - int(duration_plot_seconds // 2 * sampling_rate))
    elif first_anomaly is not None:
        start_cond = max(0, first_anomaly - pre_seconds * sampling_rate)
    else:
        print("No anomalies found.")
        return None

    samples_to_plot = int(duration_plot_seconds * sampling_rate)
    time_series_np = time_series[:, start_cond:(start_cond + samples_to_plot)]
    time = np.linspace(0, duration_plot_seconds, time_series_np.shape[1])

    fig, ax = plt.subplots(figsize=(15, 8))
    offset = np.max(np.abs(time_series_np)) * 2

    for channel_index in range(time_series_np.shape[0]):
        ax.plot(
            time,
            time_series_np[channel_index, :] + channel_index * offset,
            color="black",
            linewidth=0.5,
        )
        y_position = np.mean(time_series_np[channel_index, :] + channel_index * offset)
        ax.text(
            duration_plot_seconds + 1,
            y_position,
            f"Ch {channel_index + 1}",
            fontsize=8,
            color="black",
            verticalalignment="center",
        )

        for anomaly in anomalies[channel_index]:
            anomaly_start = anomaly[0]
            anomaly_end = anomaly_start + window_size
            anomaly_start_time = (anomaly_start - start_cond) / sampling_rate
            anomaly_end_time = (anomaly_end - start_cond) / sampling_rate
            # Only plot if at least part of the anomaly window is in the plot range
            if (anomaly_end_time >= 0) and (anomaly_start_time <= duration_plot_seconds):
                ax.axvline(x=max(0, anomaly_start_time), color='red', linestyle='--', alpha=0.7)
                ax.axvline(x=min(duration_plot_seconds, anomaly_end_time), color='blue', linestyle='--', alpha=0.7)
                ax.axvspan(
                    max(0, anomaly_start_time),
                    min(duration_plot_seconds, anomaly_end_time),
                    color='red',
                    alpha=0.1
                )

    anomalies_in_window = []
    for ch, ch_anoms in anomalies.items():
        for anomaly in ch_anoms:
            anomaly_start = anomaly[0]
            anomaly_end = anomaly_start + window_size
            if (anomaly_end >= start_cond) and (anomaly_start < start_cond + samples_to_plot):
                anomalies_in_window.append(f"Ch{ch+1}:{anomaly_start}")

    if anomalies_in_window:
        title_str = "Anomalies in window: " + ", ".join(anomalies_in_window)
    else:
        title_str = "No anomalies in plotted window"
    ax.set_title(title_str)
    legend_elements = [
        Line2D([0], [0], color='red', linestyle='--', label='Anomaly Start'),
        Line2D([0], [0], color='blue', linestyle='--', label='Anomaly End'),
        Line2D([0], [0], color='red', alpha=0.1, linewidth=10, label='Anomaly Window')
    ]
    ax.legend(handles=legend_elements, loc='upper right')
    ax.set_xlabel("Time (seconds)")
    ax.set_xlim(0, duration_plot_seconds)
    ax.set_ylabel("Amplitude + Offset")
    plt.tight_layout()
    plt.savefig(save_path)
    plt.close()

    print(f"First anomaly detected at sample: {first_anomaly}")
    return first_anomaly

start_sample = (10629 - 20) * 256
duration_plot_seconds = 90

plot_anomalies_time_series(
    signals,
    anomalies_per_channel,
    window_size=window_size,
    stride=stride,
    save_path="graphs/time_series/anomalies_time_series_10619s.png",
    duration_plot_seconds=duration_plot_seconds,
    pre_seconds=5,
    center_seconds=None,
    start_sample=start_sample
)