Neural_Processing/auto_pre_car.py at main · DanielPerodin/Neural_Processing · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
import re
from neo.io import NeuralynxIO
import pandas as pd
import numpy as np
import os
import matplotlib.pyplot as plt

data_folder = r"data"

print("Loading Neuralynx data...")
reader = NeuralynxIO(dirname=data_folder)
blk = reader.read_block(lazy=False)

print("Extracting signals...")
analogsignals = blk.segments[0].analogsignals

signal = analogsignals[0]
raw_data = np.array(signal)

names = signal.array_annotations['channel_names']
order = np.argsort([int(name.replace('CSC','')) for name in names])

raw_data = raw_data[:, order]
channel_names = names[order]

print("Corrected channel order:", channel_names)

print(f"Raw data shape: {raw_data.shape}")

raw_data = raw_data * 3.05e-8 * 1e6

sampling_rate = float(analogsignals[0].sampling_rate)
time = analogsignals[0].times.rescale('s').magnitude

print("Data loaded successfully.")

signal_rectified = np.abs(raw_data)
sliding_window = 100000
print("Computing moving average (this may take a while)...")
signal_moving_mean = (pd.DataFrame(signal_rectified).rolling(window=sliding_window, center=True, min_periods=1).mean().to_numpy())
print("Moving average complete.")
downsample_window = 1000
signal_downsample = signal_moving_mean[::downsample_window]
time_downsample = time[::downsample_window]
raw_downsampled = raw_data[::downsample_window]

plot_channels = list(range(raw_data.shape[1]))
signal_diff = np.zeros((len(signal_downsample), len(plot_channels), len(plot_channels)))

y_max_raw_downsample = np.max(raw_downsampled[:, plot_channels])

plt.figure()
plt.title("Downsampled Raw (uV)")

for i in range(len(plot_channels)):
    plt.subplot(len(plot_channels), 1, i+1)
    plt.plot(time_downsample, raw_downsampled[:, plot_channels[i]])
    plt.ylabel(f"Ch {plot_channels[i]}")
    plt.ylim(-0.1*y_max_raw_downsample, y_max_raw_downsample)

plt.savefig("downsampled_raw.png", dpi=200)
plt.close()

""" Compute pairwise channel differences
for plot_index in range(len(plot_channels)):
    for plot_diff_index in range(len(plot_channels)):
        signal_diff[:, plot_index, plot_diff_index] = np.abs(
            signal_downsample[:, plot_index] -
            signal_downsample[:, plot_diff_index]
        )

y_max_downsample = np.max(signal_downsample[:, plot_channels])
y_max_diff = np.max(signal_diff[:, plot_channels][:, :, plot_channels]) """

""" Plot rectified downsampled signals
plt.figure()
plt.title("Downsampled Rectified (uV)")
for plot_index in range(len(plot_channels)):
    plt.subplot(len(plot_channels), 1, plot_index+1)
    plt.plot(time_downsample, signal_downsample[:, plot_channels[plot_index]])
    plt.ylabel(f"Ch {plot_channels[plot_index]+1}")
    plt.ylim(-0.1*y_max_downsample, y_max_downsample)
plt.tight_layout()
plt.show() """

""" Plot pairwise difference matrix (very large plot!)
plt.figure()
for plot_index in range(len(plot_channels)):
    for plot_diff_index in range(len(plot_channels)):
        plt.subplot(
            len(plot_channels),
            len(plot_channels),
            plot_index*len(plot_channels) + plot_diff_index + 1

        plt.plot(time_downsample, signal_diff[:, plot_index, plot_diff_index])
        plt.ylabel(f"Diff Ch {plot_index+1}")
        plt.ylim(-0.1*y_max_downsample, y_max_downsample)
plt.tight_layout()
plt.show() """

signal_sorted = np.sort(signal_downsample, axis=0)[::-1]
n = round(signal_sorted.shape[0] / 100)
signal_peaks = signal_sorted[:n, :]
signal_peaks_mean = np.mean(signal_peaks, axis=0)
signal_peaks_mean_sorted_index = np.argsort(signal_peaks_mean)[::-1]
signal_peaks_mean_sorted = signal_peaks_mean[signal_peaks_mean_sorted_index]

correlation_matrix = np.corrcoef(signal_downsample, rowvar=False)
print("Correlation matrix done.")
high_correlation = correlation_matrix > 0.9
isZero = (high_correlation == 0)
[row, col] = np.nonzero(isZero)

non_group_0 = list(signal_peaks_mean_sorted_index)
group_1 = [non_group_0[0]]
non_group_1 = []

for ch in non_group_0[1:]:
    if high_correlation[non_group_0[0], ch]:
        group_1.append(ch)
    else:
        non_group_1.append(ch)

print("Group 1:", np.array(group_1) + 1)

non_group_2 = []
group_2 = []

if len(non_group_1) > 1:
    group_2 = [non_group_1[0]]
    for ch in non_group_1[1:]:
        if high_correlation[non_group_1[0], ch]:
            group_2.append(ch)
        else:
            non_group_2.append(ch)
else:
    group_2 = non_group_1

print("Group 2:", np.array(group_2) + 1)

non_group_3 = []
group_3 = []

if len(non_group_2) > 1:
    group_3 = [non_group_2[0]]
    for ch in non_group_2[1:]:
        if high_correlation[non_group_2[0], ch]:
            group_3.append(ch)
        else:
            non_group_3.append(ch)
else:
    group_3 = non_group_2

print("Group 3:", np.array(group_3) + 1)

print("\nFinal channel groupings")
print("Group 1:", np.array(group_1) + 1)
print("Group 2:", np.array(group_2) + 1)
print("Group 3:", np.array(group_3) + 1)