OSCR

DBSsync: combining intracranial and multimodal data to investigate new biomarkers in Parkinson's disease.

Code ↔ Paper

15 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 15 matches
  1. [1] § Methods › Cardiac artifact removal from Percept-LFP data ↔ functions/plotting.py, lines 225–303 · score 0.83 · band pass filtered, 0.5–60 Hz, manually override, external ECG channel, peak detection, 0.5 Hz
  2. [2] § Methods › Cardiac artifact removal from Percept-LFP data ↔ functions/ecg_cleaning.py, lines 264–323 · score 0.83 · band pass filtered, 0.5–60 Hz, peak detection, ECG channel, algorithm, override
  3. [3] § Methods › Synchronization of intracranial Percept-LFP and external recordings ↔ functions/find_artifacts.py, lines 51–100 · score 0.82 · threshold window, bipolar electrode, external channel, sharp, external recordings, crossing
  4. [4] § Methods › Cardiac artifact removal from Percept-LFP data ↔ DBSsync_main.py, lines 1009–1122 · score 0.79 · cardiac artifact removal, power spectrum, cleaned channel, ECG artifact, reconstruction, epochs
  5. [5] § Methods › Cardiac artifact removal from Percept-LFP data ↔ functions/ecg_cleaning.py, lines 1–44 · score 0.73 · QRS complex, template subtraction, fitted, Stam, singular, components
  6. [6] § Methods › Cardiac artifact removal from Percept-LFP data ↔ DBSsync_main.py, lines 1009–1122 · score 0.69 · cardiac artifact removal, ECG artifacts, fitted, QRS, subtraction, reconstructed
  7. [7] § Methods › Synchronization of intracranial Percept-LFP and external recordings ↔ functions/interactive.py, lines 129–207 · score 0.68 · pass filter, external ECG channel, cardiac artifacts, detrend, doesn, threshold
  8. [8] § Methods › Validation of cardiac artifact removal ↔ functions/ecg_cleaning.py, lines 264–323 · score 0.68 · externally recorded ECG, LFP signal, pass filtered, validate, window, artifact
  9. [9] § Methods › Main interface and compatible file formats ↔ functions/io.py, lines 1–50 · score 0.67 · MNE python, compatibility, POLY5, FIF, interface, GUI
  10. [10] § Methods › Validation of cardiac artifact removal ↔ functions/ecg_cleaning.py, lines 61–183 · score 0.67 · avoid stimulation pulses, manual override, peak detection, polarity, ECG, artifacts
  11. [11] § Methods › Cardiac artifact removal from Percept-LFP data ↔ functions/ecg_cleaning.py, lines 985–1116 · score 0.62 · power spectrum, cleaned channel, ECG artifact, epochs, overlap, template
  12. [12] § Methods › Cardiac artifact removal from Percept-LFP data ↔ functions/ecg_cleaning.py, lines 590–633 · score 0.61 · pass detection, refined, segment, epochs, correlation, threshold
  13. [13] § Results › Stimulation artifacts are reproducible across patients ↔ functions/find_artifacts.py, lines 1–26 · score 0.59 · identifiable artifacts, bipolar electrode, intracranial recordings, external recordings, stimulation, DBS
  14. [14] § Methods › Induction of DBS synchronization artifacts ↔ functions/find_artifacts.py, lines 51–100 · score 0.56 · short pulses, 0–1, ramp, external recordings, amplitude, stimulation
  15. [15] § Methods › Validation of cardiac artifact removal › Beta/theta power preservation (BPP/TPP) ↔ functions/ecg_cleaning.py, lines 985–1116 · score 0.52 · Power spectral density, Welch, PSD, overlap, window, signals

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,210 lines · 48 KB · MIT · 7 matches

  1. """
  2. This module contains functions for cleaning ECG data using various methods as
  3. described in Stam et al., 2023.
  5. There are three main methods implemented:
  6. 1. Interpolation Method: This method uses simple linear interpolation to remove artifacts
  7. from the ECG signal, i.e. linearly interpolated over the R-peaks found.
  8. 2. Template Subtraction Method: This method creates a QRS template based on the
  9. R-peaks found and subtracts it from the ECG signal, using a linear fit to adjust
  10. the template to the raw data.
  11. 3. Singular Value Decomposition (SVD) Method: This method uses SVD to decompose
  12. the raw contaminated signal into components, allowing for the removal of
  13. artifacts by selecting the most significant components.
  15. The detection of R-peaks is done using either only the LFP channel or an
  16. additional external ECG channel synchronized with the LFP channel.
  17. The R-peaks are detected with the function scipy.signal.find_peaks, with
  18. a threshold-based at 95th percentile of the signal amplitude, and a minimum
  19. distance of 500 ms between peaks. The polarity of the R-peaks is determined
  20. by comparing the mean amplitude of the detected peaks in both orientations
  21. (positive and negative). The orientation with the higher mean absolute amplitude
  22. is chosen as the orientation of the QRS complexes.
  23. There is the possibility to manually override some parameters for R-peak detection,
  24. such as the polarity, start and end times for cleaning, periods to exclude peaks from,
  25. and the detection threshold.
  27. Reference paper:
  28. Stam M.J., van Wijk B.C.M., Sharma P., Beudel M., Piña-Fuentes D.A.,
  29. de Bie R.M.A., Schuurman P.R., Neumann W.J., Buijink A.W.G. (2023)
  30. A comparison of methods to suppress electrocardiographic artifacts in local
  31. field potential recordings. Clin Neurophysiol. doi: 10.1016/j.clinph.2022.11.011
  33. This module contains the following functions:
  34. - `find_r_peaks`: Finds R-peaks in the LFP channel using either only the LFP channel itself,
  35. or an additional external ECG channel synchronized with the LFP channel.
  36. - `find_r_peaks_based_on_ext_ecg`: Finds R-peaks in the LFP channel based on an external ECG channel.
  37. - `find_r_peaks_in_lfp_channel`: Finds R-peaks in the LFP channel using the LFP channel itself.
  38. - `start_ecg_cleaning_interpolation`: Starts the ECG cleaning process using the interpolation method.
  39. - `clean_ecg_interpolation`: Cleans the ECG signal using the interpolation method.
  40. - `start_ecg_cleaning_template_sub`: Starts the ECG cleaning process using the template subtraction method.
  41. - `clean_ecg_template_sub`: Cleans the ECG signal using the template subtraction method.
  42. - `start_ecg_cleaning_svd`: Starts the ECG cleaning process using the SVD method.
  43. - `clean_ecg_svd`: Cleans the ECG signal using the SVD method.
  44. """
  46. #######################################################################
  47. # ECG CLEANING FUNCTIONS #
  48. #######################################################################
  50. from PyQt5.QtWidgets import QMessageBox, QDialog, QVBoxLayout, QLabel, QLineEdit, QComboBox, QHBoxLayout, QPushButton
  51. import numpy as np
  52. import scipy
  53. import scipy.signal
  54. import mne
  55. import ast
  57. from functions.utils import get_start_end_times, find_similar_sample
  58. from functions.classes import PlotWindow
  61. def manual_override(self):
  62. # Create dialog box for manual parameter input
  63. dialog = QDialog()
  64. dialog.setWindowTitle("Set R-peak Detection Parameters")
  66. layout = QVBoxLayout()
  68. # Create combo box for polarity selection
  69. polarity_layout = QHBoxLayout()
  70. polarity_label = QLabel("R-peak polarity (LFP):")
  71. combo_polarity = QComboBox()
  72. combo_polarity.addItems(["None", "Down", "Up"])
  73. # combo_polarity.setCurrentText("None")
  74. combo_polarity.setCurrentText(self.r_peak_polarity_lfp)
  75. polarity_layout.addWidget(polarity_label)
  76. polarity_layout.addWidget(combo_polarity)
  77. layout.addLayout(polarity_layout)
  79. # Line edits to adapt start/end cleaning times (avoid stimulation pulses)
  80. start_layout = QHBoxLayout()
  81. start_label = QLabel("Start cleaning time (s):")
  82. start_edit = QLineEdit()
  83. # start_edit.setPlaceholderText("None")
  84. start_edit.setText(str(self.start_cleaning_time))
  85. start_layout.addWidget(start_label)
  86. start_layout.addWidget(start_edit)
  87. layout.addLayout(start_layout)
  88. end_layout = QHBoxLayout()
  89. end_label = QLabel("End cleaning time (s):")
  90. end_edit = QLineEdit()
  91. # end_edit.setPlaceholderText("None")
  92. end_edit.setText(str(self.end_cleaning_time))
  93. end_layout.addWidget(end_label)
  94. end_layout.addWidget(end_edit)
  95. layout.addLayout(end_layout)
  97. # Line edit for periods to exclude peaks from (avoid strong artifacts)
  98. exclusion_layout = QHBoxLayout()
  99. exclusion_label = QLabel("Exclusion period (s) as tuples:")
  100. exclusion_edit = QLineEdit()
  101. # exclusion_edit.setPlaceholderText("None")
  102. exclusion_edit.setText(str(self.exclusion_periods))
  103. exclusion_layout.addWidget(exclusion_label)
  104. exclusion_layout.addWidget(exclusion_edit)
  105. layout.addLayout(exclusion_layout)
  107. # Combo box for R-peak detection threshold. Default to 95%
  108. threshold_layout = QHBoxLayout()
  109. threshold_label = QLabel("R-peak detection threshold (%):")
  110. combo_r_peak_threshold = QComboBox()
  111. combo_r_peak_threshold.addItems(["95", "96", "97", "98", "99"])
  112. combo_r_peak_threshold.setCurrentText(str(self.detection_threshold))
  113. threshold_layout.addWidget(threshold_label)
  114. threshold_layout.addWidget(combo_r_peak_threshold)
  115. layout.addLayout(threshold_layout)
  117. # OK / Cancel buttons
  118. button_layout = QHBoxLayout()
  119. ok_button = QPushButton("OK")
  120. cancel_button = QPushButton("Cancel")
  121. button_layout.addWidget(ok_button)
  122. button_layout.addWidget(cancel_button)
  123. layout.addLayout(button_layout)
  125. dialog.setLayout(layout)
  127. # Button connections
  128. def on_ok():
  129. try:
  130. r_peak_polarity_lfp = combo_polarity.currentText()
  131. if r_peak_polarity_lfp == "None":
  132. r_peak_polarity_lfp = None
  134. start_text = start_edit.text().strip()
  135. if start_text == "None":
  136. start_text = None
  137. start_cleaning_time = (float(start_text) if start_text is not None else None)
  139. end_text = end_edit.text().strip()
  140. if end_text == "None":
  141. end_text = None
  142. end_cleaning_time = (float(end_text) if end_text is not None else None)
  144. exclusion_text = exclusion_edit.text().strip()
  145. if exclusion_text == "None":
  146. exclusion_text = None
  147. if exclusion_text:
  148. try:
  149. exclusion_periods = ast.literal_eval(exclusion_text)
  150. # Check it's a list of tuples of floats
  151. if (isinstance(exclusion_periods, list) and
  152. all(isinstance(t, tuple) and len(t) == 2 for t in exclusion_periods)):
  153. exclusion_periods = [(float(a), float(b)) for a, b in exclusion_periods]
  154. else:
  155. raise ValueError
  156. except Exception:
  157. QMessageBox.warning(dialog, "Invalid Input",
  158. "Please enter exclusion periods in the format [(6.2, 8.8), (101.45, 123.65)].")
  159. return
  160. else:
  161. exclusion_periods = None
  163. # store as instance attributes for later use
  164. self.r_peak_polarity_lfp = r_peak_polarity_lfp
  165. self.start_cleaning_time = start_cleaning_time
  166. self.end_cleaning_time = end_cleaning_time
  167. self.exclusion_periods = exclusion_periods
  168. self.detection_threshold = int(combo_r_peak_threshold.currentText() or 95)
  170. dialog.accept()
  171. except ValueError:
  172. QMessageBox.warning(dialog, "Invalid Input", "Please enter valid numbers for times.")
  174. ok_button.clicked.connect(on_ok)
  175. cancel_button.clicked.connect(dialog.reject)
  177. # Show dialog
  178. if dialog.exec_() == QDialog.Accepted:
  179. print(f"Polarity: {self.r_peak_polarity_lfp}, "
  180. f"Start: {self.start_cleaning_time}, End: {self.end_cleaning_time}, "
  181. f"Ignoring: {self.exclusion_periods}, Threshold: {self.detection_threshold}%")
  182. else:
  183. return # User canceled
  185. ###############################################################################
  187. #######################################################################
  188. ######### FINDING R-PEAKS #########
  189. #######################################################################
  191. def find_r_peaks(self):
  192. """
  193. Find R-peaks in the LFP channel using either only the LFP channel itself,
  194. or an additional external ECG channel synchronized with the LFP channel.
  195. This function will be called when the user clicks the "Find R-peaks" button.
  196. """
  198. if self.config['NoSync'] == True:
  199. full_data = self.dataset_intra.raw_data.get_data()[
  200. self.dataset_intra.selected_channel_index_ecg
  201. ]
  203. else:
  204. full_data = self.dataset_intra.synced_data.get_data()[
  205. self.dataset_intra.selected_channel_index_ecg
  206. ]
  208. nb_points = len(full_data)
  209. times = np.arange(nb_points) * (1.0 / self.dataset_intra.sf)
  211. if self.dataset_extra.selected_channel_name_ecg is not None:
  212. # an external ECG channel was selected and will therefore be used for R-peak detection
  213. # Use external ECG channel to find R-peaks
  214. final_peaks, polarity, mean_epoch = find_r_peaks_based_on_ext_ecg(
  215. self, full_data, times, self.detection_threshold,
  216. window_artifact = [-0.5, 0.5]
  217. )
  218. # Add a message window to inform the user about the result of the peak detection:
  219. QMessageBox.information(
  220. self,
  221. "R-Peak Detection",
  222. f"{len(final_peaks)} R-peaks have been detected in the LFP channel using the ECG channel provided.",
  223. QMessageBox.Ok
  224. )
  225. else:
  226. # Use only LFP channel to find R-peaks
  227. final_peaks, polarity, mean_epoch = find_r_peaks_in_lfp_channel(
  228. self, full_data, times, self.detection_threshold,
  229. window = [-0.5, 0.5]
  230. )
  231. # Add a message window to inform the user about the result of the peak detection:
  232. QMessageBox.information(
  233. self,
  234. "R-Peak Detection",
  235. f"{len(final_peaks)} R-peaks have been detected in the LFP channel alone, using a threshold of {self.detection_threshold}%.",
  236. QMessageBox.Ok
  237. )
  239. # Check if peak detection was successful
  240. if len(final_peaks) == 0:
  241. # Disable cleaning buttons since we don't have valid peaks
  242. self.btn_start_ecg_cleaning_interpolation.setEnabled(False)
  243. self.btn_start_ecg_cleaning_template_sub.setEnabled(False)
  244. self.btn_start_ecg_cleaning_svd.setEnabled(False)
  245. QMessageBox.warning(
  246. self,
  247. "R-Peak Detection",
  248. f"R-peak detection failed, no peaks were detected",
  249. QMessageBox.Ok
  250. )
  251. return
  253. else:
  254. self.final_peaks = final_peaks
  255. self.polarity = polarity
  256. self.mean_epoch = mean_epoch
  258. self.btn_start_ecg_cleaning_interpolation.setEnabled(True)
  259. self.btn_start_ecg_cleaning_template_sub.setEnabled(True)
  260. self.btn_start_ecg_cleaning_svd.setEnabled(True)
  264. def find_r_peaks_based_on_ext_ecg(
  265. self,
  266. full_data: np.ndarray,
  267. times: np.ndarray,
  268. detection_threshold: int,
  269. window_artifact: list = [-0.5, 0.5]
  270. ):
  271. #### PREDETERMINE R-PEAKS TIMESTAMPS USING ECG CHANNEL ####
  272. """
  273. The externally recorded ECG signal is used to predetermine
  274. the timestamps of the R-peaks. The ECG signal is z-scored ((x-l)/r)
  275. over the entire recording and the function findpeaks was used to search
  276. for R-peaks with a specific height (95th percentile) and at a
  277. specific inter-peak distance (minimally 500 ms).
  278. The algorithm accounts for negative QRS complexes
  279. by repeating this procedure after multiplying the signal with
  280. -1. For both orientations of the LFP signal, the values of the
  281. peaks were averaged and the peaks with the highest mean
  282. determined the orientation of the QRS complexes.
  283. """
  284. # searches for stimulation pulses to avoid them during R-peak detection
  285. last_peak_start, first_peak_end = get_start_end_times(full_data, times)
  287. # Override with user-defined times if provided
  288. if self.start_cleaning_time is not None:
  289. last_peak_start = self.start_cleaning_time
  290. if self.end_cleaning_time is not None:
  291. first_peak_end = self.end_cleaning_time
  293. # Validate crop range and correct to default if invalid
  294. if last_peak_start >= first_peak_end:
  295. # Use a reasonable portion of the signal (skip first and last 10 seconds)
  296. last_peak_start = max(10.0, times[0] + 10.0)
  297. first_peak_end = min(times[-1] - 10.0, times[-1] - 10.0)
  299. data_extra = self.dataset_extra.synced_data.get_data()[
  300. self.dataset_extra.selected_channel_index_ecg
  301. ]
  303. # Apply 0.5 Hz-60Hz band-pass filter to ECG data. The low-pass can be changed in the config file.
  304. # b, a = scipy.signal.butter(1, 0.05, "highpass")
  305. b, a = scipy.signal.butter(1, 0.5, "highpass", fs = self.dataset_extra.sf)
  306. detrended_data = scipy.signal.filtfilt(b, a, data_extra)
  308. b2, a2 = scipy.signal.butter(
  309. N=4, # Filter order
  310. Wn=self.config["EcgLowpassFilter"],
  311. btype="lowpass",
  312. fs=self.dataset_extra.sf
  313. )
  314. ecg_data = scipy.signal.filtfilt(b2, a2, detrended_data)
  315. nb_points = len(ecg_data)
  316. timescale_extra = np.arange(nb_points) * (1.0 / self.dataset_extra.sf)
  317. # timescale_extra = np.linspace(
  318. # 0,
  319. # self.dataset_extra.synced_data.get_data().shape[1]/self.dataset_extra.sf,
  320. # self.dataset_extra.synced_data.get_data().shape[1]
  321. # )
  322. # end_time_extra = self.dataset_extra.synced_data.get_data().shape[1]/self.dataset_extra.sf
  323. # timescale_extra = np.arange(0, end_time_extra, 1/self.dataset_extra.sf)
  325. # Z-score the ECG signal
  326. ecg_z = (ecg_data - np.mean(ecg_data)) / np.std(ecg_data)
  328. # Define peak detection params in the ECG channel
  329. threshold = np.percentile(ecg_z, detection_threshold) # Convert to percentile for robustness
  330. min_distance_samples = int(0.5 * self.dataset_extra.sf) # 500 ms in samples
  332. # Detect peaks in original signal
  333. peaks_pos, _ = scipy.signal.find_peaks(
  334. ecg_z,
  335. height=threshold,
  336. distance=min_distance_samples
  337. )
  339. # Detect peaks in inverted signal
  340. peaks_neg, __ = scipy.signal.find_peaks(
  341. -ecg_z,
  342. height=threshold,
  343. distance=min_distance_samples
  344. )
  346. # Select better polarity based on the number of peaks detected
  347. if len(peaks_pos) >= len(peaks_neg):
  348. chosen_peaks = peaks_pos
  349. polarity_ecg = 'Positive'
  350. else:
  351. chosen_peaks = peaks_neg
  352. polarity_ecg = 'Negative'
  354. # Check if any peaks were detected in the external ECG
  355. if len(chosen_peaks) == 0:
  356. # Try with lower thresholds
  357. for fallback_threshold in [90, 85, 80, 75, 70]:
  358. threshold_fallback = np.percentile(ecg_z, fallback_threshold)
  360. peaks_pos_fb, _ = scipy.signal.find_peaks(
  361. ecg_z, height=threshold_fallback, distance=min_distance_samples
  362. )
  363. peaks_neg_fb, _ = scipy.signal.find_peaks(
  364. -ecg_z, height=threshold_fallback, distance=min_distance_samples
  365. )
  367. if len(peaks_pos_fb) > 0 or len(peaks_neg_fb) > 0:
  368. chosen_peaks = peaks_pos_fb if len(peaks_pos_fb) >= len(peaks_neg_fb) else peaks_neg_fb
  369. polarity_ecg = 'Positive' if len(peaks_pos_fb) >= len(peaks_neg_fb) else 'Negative'
  370. break
  372. # If still no peaks found
  373. if len(chosen_peaks) == 0:
  374. print("No peaks detected in external ECG with any threshold!")
  375. QMessageBox.warning(
  376. self,
  377. "External ECG Peak Detection Failed",
  378. "No R-peaks could be detected in the external ECG channel. Please check if:\n"
  379. "1. The correct ECG channel is selected\n"
  380. "2. The ECG signal quality is sufficient\n"
  381. "3. The synchronization between recordings is correct",
  382. QMessageBox.Ok
  383. )
  384. return np.array([]), 'Unknown', np.array([])
  387. # Convert R-peaks from ECG samples to seconds
  388. r_peak_times_sec = chosen_peaks / self.dataset_extra.sf
  390. # Convert times to LFP sample indices
  391. r_peaks_lfp_idx = np.round(r_peak_times_sec * self.dataset_intra.sf).astype(int)
  393. # Look for R-peaks in LFP channel based on predetermined timestamps:
  394. window_around_peaks = 20 # ±20 LFP samples = 80ms at 250 Hz
  395. max_peaks = []
  396. min_peaks = []
  398. for idx in r_peaks_lfp_idx:
  399. start = max(idx - window_around_peaks, 0)
  400. end = min(idx + window_around_peaks + 1, len(full_data))
  401. segment = full_data[start:end]
  403. if len(segment) > 0:
  404. max_peaks.append(np.max(segment))
  405. min_peaks.append(np.min(segment))
  407. # Calculate mean absolute values
  408. mean_abs_max = np.nanmean(np.abs(max_peaks)) if max_peaks else 0
  409. mean_abs_min = np.nanmean(np.abs(min_peaks)) if min_peaks else 0
  411. # Override polarity if user specified
  412. if self.r_peak_polarity_lfp is not None:
  413. polarity = self.r_peak_polarity_lfp
  414. if polarity == 'Up':
  415. mean_abs_max = 2
  416. mean_abs_min = 1
  417. else:
  418. mean_abs_max = 1
  419. mean_abs_min = 2
  421. # Choose the orientation with the higher mean absolute amplitude in the LFP channel
  422. lfp_peak_indices = []
  423. polarity = None
  424. if mean_abs_max >= mean_abs_min:
  425. polarity = 'Up'
  426. for idx in r_peaks_lfp_idx:
  427. start = idx - window_around_peaks
  428. end = idx + window_around_peaks + 1
  429. # Check signal boundaries
  430. if start < 0 or end > len(full_data):
  431. continue
  432. segment = full_data[start:end]
  433. if np.isnan(segment).any():
  434. continue
  435. local_max_idx = np.argmax(segment)
  436. peak_global_idx = start + local_max_idx
  437. lfp_peak_indices.append(peak_global_idx)
  438. else:
  439. polarity = 'Down'
  440. for idx in r_peaks_lfp_idx:
  441. start = idx - window_around_peaks
  442. end = idx + window_around_peaks + 1
  443. # Check signal boundaries
  444. if start < 0 or end > len(full_data):
  445. continue
  446. segment = full_data[start:end]
  447. if np.isnan(segment).any():
  448. continue
  449. local_min_idx = np.argmin(segment)
  450. peak_global_idx = start + local_min_idx
  451. lfp_peak_indices.append(peak_global_idx)
  453. # Remove peaks that are before last_peak_start and after first_peak_end:
  454. initial_lfp_count = len(lfp_peak_indices)
  455. lfp_peak_indices = [
  456. p for p in lfp_peak_indices if (
  457. p >= int(last_peak_start * self.dataset_intra.sf) and
  458. p <= int(first_peak_end * self.dataset_intra.sf)
  459. )]
  461. # Also remove peaks in exclusion periods if provided
  462. if self.exclusion_periods is not None:
  463. lfp_peak_indices = [
  464. p for p in lfp_peak_indices
  465. if not any(start <= (p / self.dataset_intra.sf) <= end for start, end in self.exclusion_periods)
  466. ]
  468. # Check if we have sufficient LFP peaks
  469. if len(lfp_peak_indices) == 0:
  470. print("No LFP peaks remain after time filtering!")
  471. QMessageBox.warning(
  472. self,
  473. "LFP Peak Detection Failed",
  474. "No valid R-peaks found in the LFP channel after applying time constraints. "
  475. "This might be due to timing issues between the external ECG and LFP recordings.",
  476. QMessageBox.Ok
  477. )
  478. return np.array([]), polarity, np.array([])
  480. # Plot detected R-peaks in external
  481. # first scale the ECG data to match the amplitude of the LFP channel for better visualization
  482. ptp_lfp = np.ptp(full_data) / 2
  483. ptp_ecg = np.ptp(ecg_data)
  484. factor = ptp_lfp/ptp_ecg
  485. ecg_data_scaled = ecg_data * factor
  487. self.canvas_detected_peaks.setEnabled(True)
  488. self.toolbar_detected_peaks.setEnabled(True)
  489. self.ax_detected_peaks.clear()
  490. self.ax_detected_peaks.set_title('Detected Peaks')
  491. self.ax_detected_peaks.plot(timescale_extra, ecg_data_scaled, label='Raw ECG', alpha=0.1)
  492. self.ax_detected_peaks.plot(timescale_extra[chosen_peaks], ecg_data_scaled[chosen_peaks], 'ro', label='Detected Peaks', alpha=0.1)
  493. self.canvas_detected_peaks.draw()
  495. # Plot detected R-peaks in intracranial
  496. self.ax_detected_peaks.plot(times, full_data, label='Raw LFP', color='black')
  497. self.ax_detected_peaks.plot(
  498. np.array(times)[lfp_peak_indices],
  499. np.array(full_data)[lfp_peak_indices],
  500. 'ro', label='LFP Peaks'
  501. )
  502. self.ax_detected_peaks.legend()
  503. self.canvas_detected_peaks.draw()
  505. # Estimate HR and display it in label
  506. peak_intervals = np.diff(lfp_peak_indices) / self.dataset_intra.sf # Convert to seconds
  507. hr = 60 / np.mean(peak_intervals) if len(peak_intervals) > 0 else 0
  508. self.label_heart_rate_lfp.setText(f'Heart rate: {hr:.1f} bpm')
  510. # Define epoch window
  511. sf_lfp = self.dataset_intra.sf
  512. pre_samples = int(abs(window_artifact[0]) * sf_lfp)
  513. post_samples = int(window_artifact[1] * sf_lfp)
  514. epoch_length = pre_samples + post_samples # Total length of each epoch
  515. # time = np.linspace(window_artifact[0], window_artifact[1], epoch_length) # Time in seconds
  516. time_epoch = np.arange(-pre_samples, post_samples) / sf_lfp # Time in seconds
  518. epochs = [] # Store extracted heartbeats
  520. for peak in lfp_peak_indices:
  521. start = peak - pre_samples
  522. end = peak + post_samples
  523. if np.isnan(full_data[start:end]).any():
  524. print(f"Skipping peak at {peak} due to NaNs in the epoch")
  525. continue
  526. if (
  527. start >= last_peak_start*self.dataset_intra.sf
  528. ) and (
  529. end < first_peak_end*self.dataset_intra.sf
  530. ): # Ensure we don't take the peaks that are in the stimulation pulses
  531. epochs.append(full_data[start:end])
  533. epochs = np.array(epochs)
  535. # Check if we have valid epochs
  536. if len(epochs) == 0:
  537. print("No valid epochs extracted from LFP peaks!")
  538. QMessageBox.warning(
  539. self,
  540. "LFP Epoch Extraction Failed",
  541. "No valid epochs could be extracted from the detected LFP peaks. "
  542. "This might be due to peaks being too close to signal boundaries.",
  543. QMessageBox.Ok
  544. )
  545. return np.array([]), polarity, np.array([])
  547. # Compute average heartbeat template
  548. mean_epoch = np.nanmean(epochs, axis=0)
  550. # Additional check for mean_epoch validity
  551. if np.isnan(mean_epoch).all() or len(mean_epoch) == 0:
  552. print("LFP mean epoch is invalid (all NaN or empty)!")
  553. QMessageBox.warning(
  554. self,
  555. "LFP Template Creation Failed",
  556. "Could not create a valid ECG template from the detected LFP peaks.",
  557. QMessageBox.Ok
  558. )
  559. return np.array([]), polarity, np.array([])
  561. # Plot the detected ECG epochs
  562. self.canvas_ecg_artifact.setEnabled(True)
  563. self.toolbar_ecg_artifact.setEnabled(True)
  564. self.ax_ecg_artifact.clear()
  565. self.ax_ecg_artifact.set_title("Detected ECG epochs")
  567. for epoch in epochs:
  568. self.ax_ecg_artifact.plot(time_epoch, epoch, color='gray', alpha=0.3)
  570. self.ax_ecg_artifact.plot(
  571. time_epoch,
  572. mean_epoch,
  573. color='black',
  574. linewidth=2,
  575. label='Average ECG Template'
  576. )
  577. self.ax_ecg_artifact.set_xlabel("Time (s)")
  578. self.ax_ecg_artifact.set_ylabel("Amplitude")
  579. self.ax_ecg_artifact.legend()
  580. self.canvas_ecg_artifact.draw()
  582. ## store start and end times of cleaning
  583. self.after_first_stim_pulses = last_peak_start
  584. self.before_last_stim_pulses = first_peak_end
  586. return lfp_peak_indices, polarity, mean_epoch
  590. def find_r_peaks_in_lfp_channel(
  591. self,
  592. full_data: np.ndarray,
  593. times: np.ndarray,
  594. detection_threshold: int,
  595. window = [-0.5, 0.5]
  596. ):
  597. """
  598. The LFP signal itself is used to find the R-peaks. It is cropped in 1s
  599. segments and the function findpeaks was used to search for R-peaks in each segment.
  600. This is repeated with signal multiplied by -1 to account for negative QRS complexes.
  601. The polarity for which more peaks were detected is chosen as the orientation
  602. of the QRS complexes in LFP channel and a second-pass detection is performed:
  603. epochs are created around each detected R-peak (-0.5 - +0.5s around) and averaged
  604. to generate a first QRS template.
  605. Template correlation is used to refine the R-peak locations by searching
  606. for the local maxima within each epoch that has the highest correlation
  607. with the mean QRS template. These local maxima are used to create a better
  608. QRS template by averaging epochs around these new R-peak locations.
  609. Template correlation is applied a second time using this improved QRS template
  610. to finalize R-peak locations.
  611. """
  612. sf_lfp = round(self.dataset_intra.sf)
  614. # Define period to look for peaks (avoid stimulation pulses if present):
  615. if self.config['NoSync'] == True:
  616. start_idx = 0
  617. end_idx = len(full_data)
  618. # Override with user-defined times if provided
  619. if self.start_cleaning_time is not None:
  620. start_idx = int(self.start_cleaning_time * sf_lfp)
  621. if self.end_cleaning_time is not None:
  622. end_idx = int(self.end_cleaning_time * sf_lfp)
  624. last_peak_start = start_idx / sf_lfp
  625. first_peak_end = end_idx / sf_lfp
  626. else:
  627. last_peak_start, first_peak_end = get_start_end_times(full_data, times)
  629. # Override with user-defined times if provided
  630. if self.start_cleaning_time is not None:
  631. last_peak_start = self.start_cleaning_time
  632. if self.end_cleaning_time is not None:
  633. first_peak_end = self.end_cleaning_time
  635. # Validate crop range
  636. if last_peak_start >= first_peak_end:
  637. # Use a reasonable portion of the signal (skip first and last 10 seconds)
  638. last_peak_start = max(10.0, times[0] + 10.0)
  639. first_peak_end = min(times[-1] - 10.0, times[-1] - 10.0)
  641. # Calculate crop indices and validate
  642. start_idx = int(last_peak_start * sf_lfp)
  643. end_idx = int(first_peak_end * sf_lfp)
  645. # Ensure indices are within bounds
  646. start_idx = max(0, min(start_idx, len(full_data) - 1))
  647. end_idx = max(start_idx + 1, min(end_idx, len(full_data)))
  649. cropped_data = full_data[start_idx:end_idx]
  651. # Check if cropped data is valid
  652. if len(cropped_data) == 0:
  653. cropped_data = full_data
  654. last_peak_start = times[0]
  655. first_peak_end = times[-1]
  657. #ecg = {'proc': {}}
  658. ns = len(cropped_data) # Number of samples in the cropped data
  660. # Segment the signal into overlapping windows
  661. dwindow = int(round(sf_lfp)) # 1s window
  662. dmove = sf_lfp # 1s step
  663. n_segments = (ns - dwindow) // dmove + 1
  665. detected_peaks_positive = [] # Store peak indices in the original timescale of the cropped_data
  666. x = np.array(
  667. [
  668. cropped_data[
  669. i * dmove: i * dmove + dwindow
  670. ] for i in range(n_segments) if i * dmove + dwindow <= ns
  671. ])
  672. # Loop through each segment and find peaks
  673. for i in range(n_segments):
  674. segment = x[i]
  675. # Skip segment if it contains any NaNs
  676. if np.isnan(segment).any():
  677. continue
  679. for threshold_pct in [80, 70, 60, 50]:
  680. peaks, _ = scipy.signal.find_peaks(
  681. segment, height=np.percentile(segment, threshold_pct), distance=sf_lfp//3
  682. )
  683. if len(peaks) > 0:
  684. break
  686. real_peaks = peaks + (i * dmove) # Convert to original timescale
  687. detected_peaks_positive.extend(real_peaks)
  689. # Repeat with reverted signal to find negative peaks
  690. detected_peaks_negative = [] # Store peak indices in the original timescale of the cropped_data
  691. x_neg = np.array(
  692. [
  693. -cropped_data[
  694. i * dmove: i * dmove + dwindow
  695. ] for i in range(n_segments) if i * dmove + dwindow <= ns
  696. ])
  697. for i in range(n_segments):
  698. segment = x_neg[i]
  699. # 6. Skip segment if it contains any NaNs
  700. if np.isnan(segment).any():
  701. continue
  703. # Try different thresholds if 90th percentile fails
  704. for threshold_pct in [80, 70, 60, 50]:
  705. peaks, _ = scipy.signal.find_peaks(
  706. segment, height=np.percentile(segment, threshold_pct), distance=sf_lfp//3
  707. )
  708. if len(peaks) > 0:
  709. break
  711. real_peaks = peaks + (i * dmove) # Convert to original timescale
  712. detected_peaks_negative.extend(real_peaks)
  714. # Find which set of peaks has more elements
  715. if len(detected_peaks_positive) >= len(detected_peaks_negative):
  716. detected_peaks = detected_peaks_positive
  717. polarity = 'Up'
  718. else:
  719. detected_peaks = detected_peaks_negative
  720. polarity = 'Down'
  722. # Override polarity if user specified
  723. if self.r_peak_polarity_lfp is not None:
  724. polarity = self.r_peak_polarity_lfp
  725. print(f"Overriding detected polarity to user-specified: {polarity}")
  726. if polarity == 'Up':
  727. detected_peaks = detected_peaks_positive
  728. else:
  729. detected_peaks = detected_peaks_negative
  731. detected_peaks = np.array(detected_peaks)
  733. # Check if any peaks were detected
  734. if len(detected_peaks) == 0:
  735. # Try the simpler fallback method
  736. detected_peaks, polarity = simple_peak_detection_fallback(cropped_data, sf_lfp)
  738. if len(detected_peaks) == 0:
  739. print("No peaks detected with any method!")
  740. QMessageBox.warning(
  741. self,
  742. "Peak Detection Failed",
  743. "No R-peaks could be detected in the signal with any method. Please check if:\n"
  744. "1. The signal contains visible ECG artifacts\n"
  745. "2. The channel selection is correct\n"
  746. "3. The signal amplitude is sufficient\n"
  747. "4. Try adjusting the threshold percentage",
  748. QMessageBox.Ok
  749. )
  750. # Return empty arrays to prevent crashes
  751. return np.array([]), 'Unknown', np.array([])
  753. # Define epoch window
  754. pre_samples = int(abs(window[0]) * sf_lfp)
  755. post_samples = int(window[1] * sf_lfp)
  756. #epoch_length = pre_samples + post_samples # Total length of each epoch
  758. epochs = [] # Store extracted heartbeats
  759. for peak in detected_peaks:
  760. start = peak - pre_samples
  761. end = peak + post_samples
  763. if start >= 0 and end < ns: # Ensure we don't go out of bounds
  764. epoch = cropped_data[start:end]
  765. if np.isnan(epoch).any():
  766. continue
  767. else:
  768. epochs.append(epoch)
  770. epochs = np.array(epochs)
  772. # Check if we have valid epochs
  773. if len(epochs) == 0:
  774. print("No valid epochs extracted from detected peaks!")
  775. QMessageBox.warning(
  776. self,
  777. "Epoch Extraction Failed",
  778. "No valid epochs could be extracted from the detected peaks. This might be due to peaks being too close to signal boundaries.",
  779. QMessageBox.Ok
  780. )
  781. # Return empty arrays to prevent crashes
  782. return np.array([]), polarity, np.array([])
  784. # Compute average heartbeat template
  785. mean_epoch = np.nanmean(epochs, axis=0)
  787. # Additional check for mean_epoch validity
  788. if np.isnan(mean_epoch).all() or len(mean_epoch) == 0:
  789. print("Mean epoch is invalid (all NaN or empty)!")
  790. QMessageBox.warning(
  791. self,
  792. "Template Creation Failed",
  793. "Could not create a valid ECG template from the detected peaks.",
  794. QMessageBox.Ok
  795. )
  796. return np.array([]), polarity, np.array([])
  798. # Temporal correlation for ECG detection
  799. # adapt in case NaNs are present:
  800. if np.isnan(cropped_data).any():
  801. cropped_data_clean = np.nan_to_num(cropped_data, nan=0.0)
  802. r = np.correlate(cropped_data_clean, mean_epoch, mode='same')
  803. else:
  804. r = np.correlate(cropped_data, mean_epoch, mode='same')
  805. threshold = np.percentile(r, 95)
  806. detected_peaks, _ = scipy.signal.find_peaks(
  807. r, height=threshold, distance=sf_lfp//2
  808. )
  810. # Second pass for refining detection
  811. refined_template = np.nanmean(
  812. [cropped_data[
  813. p - dwindow//2 : p + dwindow//2
  814. ] for p in detected_peaks if p - dwindow//2 > 0 and p + dwindow//2 < ns
  815. ], axis=0
  816. )
  818. if np.isnan(cropped_data).any():
  819. cropped_data_clean = np.nan_to_num(cropped_data, nan=0.0)
  820. r2 = np.correlate(cropped_data_clean, refined_template, mode='same')
  821. else:
  822. r2 = np.correlate(cropped_data, refined_template, mode='same')
  823. threshold2 = np.percentile(r2, detection_threshold)
  824. final_peaks, _ = scipy.signal.find_peaks(
  825. r2, height=threshold2, distance=sf_lfp//2
  826. )
  828. # Adjust the final peaks to the original data scale
  829. final_peaks = final_peaks + start_idx
  831. # Remove peaks in exclusion periods if provided
  832. if self.exclusion_periods is not None:
  833. final_peaks = [
  834. p for p in final_peaks
  835. if not any(start <= (p / self.dataset_intra.sf) <= end for start, end in self.exclusion_periods)
  836. ]
  838. # check that no peak is close to NaN values, if yes, remove them
  839. peaks_to_remove = []
  840. for peak in final_peaks:
  841. start = peak - pre_samples
  842. end = peak + post_samples
  843. if start < 0 or end >= len(full_data):
  844. continue
  845. if np.isnan(full_data[start:end]).any():
  846. peaks_to_remove.append(peak)
  848. final_peaks = [p for p in final_peaks if p not in peaks_to_remove]
  850. # plot the detected peaks
  851. self.canvas_detected_peaks.setEnabled(True)
  852. self.toolbar_detected_peaks.setEnabled(True)
  853. self.ax_detected_peaks.clear()
  854. self.ax_detected_peaks.set_title('Detected Peaks')
  855. self.ax_detected_peaks.plot(times, full_data, label='Raw LFP', color='black')
  856. self.ax_detected_peaks.plot(
  857. np.array(times)[final_peaks],
  858. np.array(full_data)[final_peaks],
  859. 'ro', label='Detected Peaks'
  860. )
  861. self.canvas_detected_peaks.draw()
  863. # Estimate HR
  864. peak_intervals = np.diff(final_peaks) / sf_lfp # Convert to seconds
  865. hr = 60 / np.mean(peak_intervals) if len(peak_intervals) > 0 else 0
  866. self.label_heart_rate_lfp.setText(f'Heart rate: {hr} bpm')
  868. ## store start and end times of cleaning
  869. self.after_first_stim_pulses = last_peak_start
  870. self.before_last_stim_pulses = first_peak_end
  872. return final_peaks, polarity, mean_epoch
  876. #######################################################################
  877. ######### INTERPOLATION METHOD #########
  878. #######################################################################
  879. def start_ecg_cleaning_interpolation(self):
  880. self.ax_ecg_clean.clear()
  881. self.ax_ecg_artifact.clear()
  882. self.ax_psd.clear()
  883. """Start the ECG cleaning process using the interpolation method from Perceive toolbox."""
  884. try:
  885. clean_ecg_interpolation(self)
  887. except Exception as e:
  888. QMessageBox.critical(self, "Error", f"Failed to clean ECG: {e}")
  890. def clean_ecg_interpolation(self):
  891. if self.config['NoSync'] == True:
  892. full_data = self.dataset_intra.raw_data.get_data()[
  893. self.dataset_intra.selected_channel_index_ecg
  894. ]
  896. else:
  897. full_data = self.dataset_intra.synced_data.get_data()[
  898. self.dataset_intra.selected_channel_index_ecg
  899. ]
  901. nb_points = len(full_data)
  902. times = np.arange(nb_points) * (1.0 / self.dataset_intra.sf)
  904. ############################################################################
  905. # prepare a copy of the full data to store the cleaned data
  906. clean_data = np.copy(full_data)
  907. ns = len(full_data)
  909. #### INTERPOLATE DATA AT EACH R-PEAK FOUND ####
  910. # Remove artifacts (simple interpolation)
  911. for p in self.final_peaks:
  912. clean_data[max(0, p - 5): min(ns, p + 5)] = np.nan # NaN out artifacts
  913. clean_data = np.interp(
  914. np.arange(ns), np.arange(ns)[~np.isnan(clean_data)],
  915. clean_data[~np.isnan(clean_data)]
  916. )
  918. if self.dataset_intra.selected_channel_index_ecg == 0:
  919. self.dataset_intra.cleaned_ecg_left = clean_data
  920. elif self.dataset_intra.selected_channel_index_ecg == 1:
  921. self.dataset_intra.cleaned_ecg_right = clean_data
  923. # plot an overlap of the raw and cleaned data
  924. self.canvas_ecg_clean.setEnabled(True)
  925. self.toolbar_ecg_clean.setEnabled(True)
  926. self.ax_ecg_clean.clear()
  927. self.ax_ecg_clean.set_title("Cleaned ECG Signal")
  928. self.ax_ecg_clean.plot(times,full_data, label='Raw data')
  929. self.ax_ecg_clean.plot(times,clean_data, label='Cleaned data')
  930. self.ax_ecg_clean.set_xlabel("Time (s)")
  931. self.ax_ecg_clean.set_ylabel("Amplitude")
  932. self.ax_ecg_clean.legend()
  933. self.canvas_ecg_clean.draw()
  935. # Plot an overlap of the power spectrum using welch's method:
  936. n_fft = int(round(self.dataset_intra.sf))
  937. n_overlap=int(round(self.dataset_intra.sf)/2)
  939. # make sure that stimulation pulses are not included in the PSD calculation
  940. start_index = int(self.after_first_stim_pulses * self.dataset_intra.sf)
  941. end_index = int(self.before_last_stim_pulses * self.dataset_intra.sf)
  943. psd_raw, freqs_raw = mne.time_frequency.psd_array_welch(
  944. full_data[start_index: end_index],self.dataset_intra.sf,fmin=0,
  945. fmax=125,n_fft=n_fft,
  946. n_overlap=n_overlap)
  948. psd_clean, freqs_clean = mne.time_frequency.psd_array_welch(
  949. clean_data[start_index: end_index],self.dataset_intra.sf,fmin=0,
  950. fmax=125,n_fft=n_fft,
  951. n_overlap=n_overlap)
  953. self.canvas_psd.setEnabled(True)
  954. self.toolbar_psd.setEnabled(True)
  955. self.ax_psd.clear()
  956. self.ax_psd.plot(
  957. freqs_raw, np.log(psd_raw), color='blue', label='PSD raw channel'
  958. )
  959. self.ax_psd.plot(
  960. freqs_clean, np.log(psd_clean), color = 'orange',
  961. label='PSD cleaned channel'
  962. )
  963. self.ax_psd.legend()
  964. self.canvas_psd.draw()
  966. self.btn_confirm_cleaning.setEnabled(True) # Enable the button after cleaning
  970. #######################################################################
  971. ######### TEMPLATE SUBSTRACTION METHOD #########
  972. #######################################################################
  973. def start_ecg_cleaning_template_sub(self):
  974. self.ax_ecg_clean.clear()
  975. self.ax_ecg_artifact.clear()
  976. self.ax_psd.clear()
  977. """Start the ECG cleaning process using the template substraction method."""
  978. try:
  979. clean_ecg_template_sub(self)
  981. except Exception as e:
  982. QMessageBox.critical(self, "Error", f"Failed to clean ECG: {e}")
  985. def clean_ecg_template_sub(self):
  986. if self.config['NoSync'] == True:
  987. full_data = self.dataset_intra.raw_data.get_data()[
  988. self.dataset_intra.selected_channel_index_ecg
  989. ]
  991. else:
  992. full_data = self.dataset_intra.synced_data.get_data()[
  993. self.dataset_intra.selected_channel_index_ecg
  994. ]
  996. nb_points = len(full_data)
  997. times = np.arange(nb_points) * (1.0 / self.dataset_intra.sf)
  998. window = [-0.2, 0.2] # QRS complex window
  1000. ############################################################################
  1001. clean_data = np.copy(full_data)
  1002. ns = len(full_data)
  1004. # Create a QRS template #
  1005. pre_samples = int(abs(window[0]) * self.dataset_intra.sf)
  1006. post_samples = int(window[1] * self.dataset_intra.sf)
  1007. epoch_length = pre_samples + post_samples # Total length of each epoch
  1008. # timescale_epoch = np.linspace(
  1009. # window[0], window[1], epoch_length
  1010. # ) # Time in seconds
  1011. timescale_epoch = np.arange(-pre_samples, post_samples) / self.dataset_intra.sf # Time in seconds
  1013. epochs = [] # Store extracted heartbeats
  1015. for peak in self.final_peaks:
  1016. start = peak - pre_samples
  1017. end = peak + post_samples
  1018. epochs.append(full_data[start:end])
  1020. epochs = np.array(epochs)
  1022. # Compute average QRS template
  1023. mean_epoch = np.nanmean(epochs, axis=0)
  1024. self.canvas_ecg_artifact.setEnabled(True)
  1025. self.toolbar_ecg_artifact.setEnabled(True)
  1026. self.ax_ecg_artifact.clear()
  1027. self.ax_ecg_artifact.set_title("Detected QRS epochs")
  1029. for epoch in epochs:
  1030. self.ax_ecg_artifact.plot(
  1031. timescale_epoch, epoch, color='gray', alpha=0.3
  1032. )
  1034. self.ax_ecg_artifact.plot(
  1035. timescale_epoch, mean_epoch, color='black', linewidth=2,
  1036. label='Average QRS Template'
  1037. )
  1038. self.ax_ecg_artifact.set_xlabel("Time (s)")
  1039. self.ax_ecg_artifact.set_ylabel("Amplitude")
  1040. self.ax_ecg_artifact.legend()
  1041. self.canvas_ecg_artifact.draw()
  1043. ####################################################################
  1044. pre_samples = int(abs(window[0]) * self.dataset_intra.sf)
  1045. post_samples = int(window[1] * self.dataset_intra.sf)
  1047. for _, peak in enumerate(self.final_peaks):
  1048. raw_epoch = full_data[(peak - pre_samples):(peak + post_samples)]
  1049. # Prepare design matrix for linear fit (scale + offset)
  1050. X_template = np.vstack([mean_epoch, np.ones_like(mean_epoch)]).T
  1051. # Solve for optimal scale (a) and offset (b) using least squares
  1052. coeffs, _, _, _ = np.linalg.lstsq(X_template, raw_epoch, rcond=None)
  1053. a, b = coeffs
  1055. # Build fitted template
  1056. fitted_template = a * mean_epoch + b
  1058. # Equalize tails
  1059. complex_qrs_template, start_idx, end_idx = find_similar_sample(
  1060. fitted_template, tails=30
  1061. )
  1062. start = (peak - pre_samples) + start_idx
  1063. end = (peak - pre_samples) + end_idx
  1064. raw_epoch = full_data[start:end]
  1065. assert len(raw_epoch) == len(complex_qrs_template), "Raw epoch length does not match complex QRS template length"
  1066. clean_data[start:end] -= complex_qrs_template
  1068. if self.dataset_intra.selected_channel_index_ecg == 0:
  1069. self.dataset_intra.cleaned_ecg_left = clean_data
  1070. elif self.dataset_intra.selected_channel_index_ecg == 1:
  1071. self.dataset_intra.cleaned_ecg_right = clean_data
  1073. # plot an overlap of the raw and cleaned data
  1074. self.canvas_ecg_clean.setEnabled(True)
  1075. self.toolbar_ecg_clean.setEnabled(True)
  1076. self.ax_ecg_clean.clear()
  1077. self.ax_ecg_clean.set_title("Cleaned ECG Signal")
  1078. self.ax_ecg_clean.plot(times, full_data, label='Raw data')
  1079. self.ax_ecg_clean.plot(times, clean_data, label='Cleaned data')
  1080. self.ax_ecg_clean.set_xlabel("Time (s)")
  1081. self.ax_ecg_clean.set_ylabel("Amplitude")
  1082. self.ax_ecg_clean.legend()
  1083. self.canvas_ecg_clean.draw()
  1085. # Plot an overlap of the power spectrum using welch's method:
  1086. n_fft = int(round(self.dataset_intra.sf))
  1087. n_overlap=int(round(self.dataset_intra.sf)/2)
  1089. # make sure that stimulation pulses are not included in the PSD calculation
  1090. start_index = int(self.after_first_stim_pulses * self.dataset_intra.sf)
  1091. end_index = int(self.before_last_stim_pulses * self.dataset_intra.sf)
  1093. psd_raw, freqs_raw = mne.time_frequency.psd_array_welch(
  1094. full_data[start_index: end_index],self.dataset_intra.sf,fmin=0,
  1095. fmax=125,n_fft=n_fft,
  1096. n_overlap=n_overlap)
  1098. psd_clean, freqs_clean = mne.time_frequency.psd_array_welch(
  1099. clean_data[start_index: end_index],self.dataset_intra.sf,fmin=0,
  1100. fmax=125,n_fft=n_fft,
  1101. n_overlap=n_overlap)
  1103. self.canvas_psd.setEnabled(True)
  1104. self.toolbar_psd.setEnabled(True)
  1105. self.ax_psd.clear()
  1106. self.ax_psd.plot(
  1107. freqs_raw, np.log(psd_raw), color='blue', label='PSD raw channel'
  1108. )
  1109. self.ax_psd.plot(
  1110. freqs_clean, np.log(psd_clean), color = 'orange',
  1111. label='PSD cleaned channel'
  1112. )
  1113. self.ax_psd.legend()
  1114. self.canvas_psd.draw()
  1116. self.btn_confirm_cleaning.setEnabled(True) # Enable the button after cleaning
  1119. #######################################################################
  1120. ######### SINGULAR VALUE DECOMPOSITION METHOD #########
  1121. #######################################################################
  1123. def start_ecg_cleaning_svd(self):
  1124. """Start the ECG cleaning process using Singular Value Decomposition method."""
  1125. self.ax_ecg_clean.clear()
  1126. self.ax_ecg_artifact.clear()
  1127. self.ax_psd.clear()
  1128. try:
  1129. clean_ecg_svd(self)
  1131. except Exception as e:
  1132. QMessageBox.critical(self, "Error", f"Failed to clean ECG: {e}")
  1135. def clean_ecg_svd(self):
  1136. """
  1137. This function cleans the ECG signal using Singular Value Decomposition (SVD).
  1138. It extracts QRS templates around the R-peaks, performs SVD on these epochs,
  1139. and then reconstructs the signal using the first few singular values.
  1140. The function opens a secondary plotting window to visualize the SVD results,
  1141. so that the user can choose which components to keep to reconstruct the signal.
  1142. """
  1143. if self.config['NoSync'] == True:
  1144. self.full_data = self.dataset_intra.raw_data.get_data()[
  1145. self.dataset_intra.selected_channel_index_ecg
  1146. ]
  1147. else:
  1148. self.full_data = self.dataset_intra.synced_data.get_data()[
  1149. self.dataset_intra.selected_channel_index_ecg
  1150. ]
  1151. self.window = [-0.2, 0.2] # add an option to choose QRS or PQRST window??
  1153. # Create a QRS template
  1154. pre_samples = int(abs(self.window[0]) * self.dataset_intra.sf)
  1155. post_samples = int(self.window[1] * self.dataset_intra.sf)
  1156. self.epoch_length = pre_samples + post_samples # Total length of each epoch
  1158. epochs = [] # Store extracted heartbeats
  1160. for peak in self.final_peaks:
  1161. start = peak - pre_samples
  1162. end = peak + post_samples
  1163. epochs.append(self.full_data[start:end])
  1164. epochs = np.array(epochs) # shape: (n timepoints, n epochs)
  1166. ######### SINGULAR VALUE DECOMPOSITION ################
  1167. X = epochs.T # shape: (n epochs, n timepoints)
  1168. self.U, self.S, self.Vh = np.linalg.svd(X, full_matrices=False)
  1170. # Open the secondary plotting window for the SVD template
  1171. self.plot_window = PlotWindow(
  1172. self.process_value_from_plot, self.U, self.S, self.window,
  1173. self.epoch_length
  1174. )
  1175. self.plot_window.show()
  1180. def simple_peak_detection_fallback(cropped_data, sf_lfp):
  1181. """
  1182. Fallback method for peak detection using a simpler approach in case detection
  1183. with usual thresholds fails.
  1184. This method applies peak detection directly to the entire signal.
  1185. """
  1186. # Try peak detection on the entire signal with different parameters
  1187. for height_pct in [90, 80, 70, 60, 50, 40]:
  1188. for distance_factor in [2, 3, 4, 5]: # Different distance constraints
  1189. distance = sf_lfp // distance_factor
  1191. # Try positive peaks
  1192. height_pos = np.percentile(cropped_data, height_pct)
  1193. peaks_pos, _ = scipy.signal.find_peaks(
  1194. cropped_data, height=height_pos, distance=distance
  1195. )
  1197. # Try negative peaks
  1198. height_neg = np.percentile(-cropped_data, height_pct)
  1199. peaks_neg, _ = scipy.signal.find_peaks(
  1200. -cropped_data, height=height_neg, distance=distance
  1201. )
  1203. # Return the first successful detection with reasonable number of peaks
  1204. if len(peaks_pos) >= 3: # At least 3 peaks for meaningful analysis
  1205. return peaks_pos, 'Up'
  1206. elif len(peaks_neg) >= 3:
  1207. return peaks_neg, 'Down'
  1209. return np.array([]), 'Unknown'

ecg_cleaning.py at commit d91e8d1, under MIT · at the source

Overview

Authors: Juliette Vivien1,2, Charlotte E Stensholt1,2, Lucie Hortmann1, Merle Hendel1, Arian Memarpouri1, Roxanne Lofredi1,3, Jeroen G V Habets1,3, Alessia Cavallo1, Lucia K Feldmann1, Andrea A Kühn1,2,4,5,6
ORCID iDs: Juliette Vivien, Lucia K Feldmann, Andrea A Kühn
  1. Department of Neurology, Charité-Universitätsmedizin Berlin, Berlin, Germany
  2. Humboldt-Universität zu Berlin, Berlin School of Mind and Brain, Berlin, Germany
  3. Berlin Institute of Health (BIH), Berlin, Germany
  4. Bernstein Center for Computational Neuroscience, Humboldt-Universität, Berlin, Germany
  5. NeuroCure, Exzellenzcluster, Charité-Universitätsmedizin Berlin, Berlin, Germany
  6. DZNE, German Center for Neurodegenerative Diseases, Berlin, Germany
Institutions: Humboldt-Universität zu Berlin (Germany); Charité - Universitätsmedizin Berlin (Germany); Berlin Institute of Health at Charité - Universitätsmedizin Berlin (Germany); German Center for Neurodegenerative Diseases (Germany); Bernstein Center for Computational Neuroscience Berlin (Germany); NeuroCure Cluster of Excellence (Germany)
Journal: NPJ Parkinson's disease, volume 12, issue 1, article 151
Dates: received 28 November 2025; accepted 10 June 2026; published online 19 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41531-026-01439-z · PMID 42315529 · PMCID PMC13279824 · OpenAlex W4417115847
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.21203/rs.3.rs-8228751/v1
Status: code verified
Categories: Parkinson's (population), clinical / translational (subfield)
Methods: Spectral & time-frequency, Statistics, Machine learning, Preprocessing, Connectivity, Evoked potentials
Keywords: Biological techniques, Biomarkers, Engineering, Neurology, Neuroscience
Topic: Neurological disorders and treatments (Neurology, Medicine), according to OpenAlex
Funding: Deutscher Akademischer Austauschdienst; Deutsche Forschungsgemeinschaft (Project-ID 424778381 - TRR 295, Project-ID 424778381 – TRR 295); Berlin Institute of Health; Lundbeck Foundation (Grant Nr. R336-2020-1035)
Citations: not cited yet (Europe PMC); 17 references in the paper

Abstract

Implanted deep brain stimulation devices are now capable of chronically recording activity from intracranial brain areas during stimulation. This new type of data has the potential to increase our understanding of disease-related brain activity and its modulation in response to therapy or other types of stimuli. With the innovative approach of adaptive deep brain stimulation now clinically available, multimodal characterization of neural biomarkers becomes of utmost importance to define optimal feedback signals for adaptive brain stimulation and allow for better fine-tuning of stimulation parameters. To investigate these biomarkers, we developed DBSsync, a paradigm and an open-source Python toolbox with its graphical user interface for temporally precise synchronization of intracranial recordings with external data, allowing for multimodal research protocols. DBSsync achieves a temporal precision of 8 ms and incorporates cardiac artifact removal methods to facilitate intracranial data preprocessing, thus enabling the integration and precise synchronization of multiple brain signals with external sensors and various behavioral timeline data.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 15 matches between paragraphs and lines of code.

juliettevivien/DBSsync

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d91e8d186a5eda8a5f57897105f2487a543d7af0, 11 June 2026
Languages: Python (23)
Size: 51 files, 23 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (environment.yml, environment_linux.yml, requirements.txt)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (11 files), MNE-Python (8 files), SciPy (8 files), Matplotlib (6 files), pandas (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
25 files
At the source: github.com/juliettevivien/DBSsync

neuromodulation/perceive

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 23b12f2da3093f623cb192157e4dfc65b19a0c70, 10 September 2026
Languages: MATLAB (83), Shell (2)
Size: 150 files, 85 scripts
Software Heritage: not archived
Found in: the text, “Main interface and compatible file formats”
Holds: README, license file, tests, documentation
Not found: CITATION.cff, environment file, continuous integration
Tools: Statistics and Machine Learning Toolbox (8 files), Signal Processing Toolbox (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
87 files
At the source: github.com/neuromodulation/perceive/

Code availability

The code of the toolbox is publicly available in a Github repository and can be accessed via this link: https://github.com/juliettevivien/DBSsync.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 108 scripts, each with its path and the digest of its content;
  • 15 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Data are available conditionally through data-sharing agreements in accordance with data privacy statements signed by the patients within the legal framework of the General Data Protection Regulation of the European Union, within a time frame of 6 months. Requests should be directed to AAK or the Open Data officer.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 5 keywords, 4 funders, 16 references.

Cite

This paper

Vivien, J., Stensholt, C. E., Hortmann, L., Hendel, M., Memarpouri, A., Lofredi, R., Habets, J. G. V., Cavallo, A., Feldmann, L. K., & Kühn, A. A. (2026). DBSsync: combining intracranial and multimodal data to investigate new biomarkers in Parkinson's disease. NPJ Parkinson's disease, 12(1), 151. https://doi.org/10.1038/s41531-026-01439-z

BibTeX

@article{vivien2026dbssync,
author = {Vivien, Juliette and Stensholt, Charlotte E and Hortmann, Lucie and Hendel, Merle and Memarpouri, Arian and Lofredi, Roxanne and Habets, Jeroen G V and Cavallo, Alessia and Feldmann, Lucia K and Kühn, Andrea A},
title = {{DBSsync: combining intracranial and multimodal data to investigate new biomarkers in Parkinson's disease}},
journal = {NPJ Parkinson's disease},
year = {2026},
month = jun,
volume = {12},
number = {1},
pages = {151},
publisher = {Nature Publishing Group},
issn = {2373-8057},
doi = {10.1038/s41531-026-01439-z},
url = {https://doi.org/10.1038/s41531-026-01439-z},
pmid = {42315529},
pmcid = {PMC13279824}
}

RIS

TY - JOUR
AU - Vivien, Juliette
AU - Stensholt, Charlotte E
AU - Hortmann, Lucie
AU - Hendel, Merle
AU - Memarpouri, Arian
AU - Lofredi, Roxanne
AU - Habets, Jeroen G V
AU - Cavallo, Alessia
AU - Feldmann, Lucia K
AU - Kühn, Andrea A
TI - DBSsync: combining intracranial and multimodal data to investigate new biomarkers in Parkinson's disease
T2 - NPJ Parkinson's disease
J2 - NPJ Parkinsons Dis
PY - 2026
DA - 2026/06/19
VL - 12
IS - 1
SP - 151
SN - 2373-8057
PB - Nature Publishing Group
DO - 10.1038/s41531-026-01439-z
UR - https://doi.org/10.1038/s41531-026-01439-z
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41531-026-01439-z",
"type": "article-journal",
"title": "DBSsync: combining intracranial and multimodal data to investigate new biomarkers in Parkinson's disease",
"container-title": "NPJ Parkinson's disease",
"author": [
{
"family": "Vivien",
"given": "Juliette"
},
{
"family": "Stensholt",
"given": "Charlotte E"
},
{
"family": "Hortmann",
"given": "Lucie"
},
{
"family": "Hendel",
"given": "Merle"
},
{
"family": "Memarpouri",
"given": "Arian"
},
{
"family": "Lofredi",
"given": "Roxanne"
},
{
"family": "Habets",
"given": "Jeroen G V"
},
{
"family": "Cavallo",
"given": "Alessia"
},
{
"family": "Feldmann",
"given": "Lucia K"
},
{
"family": "Kühn",
"given": "Andrea A"
}
],
"container-title-short": "NPJ Parkinsons Dis",
"volume": "12",
"issue": "1",
"page": "151",
"DOI": "10.1038/s41531-026-01439-z",
"PMID": "42315529",
"PMCID": "PMC13279824",
"ISSN": "2373-8057",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41531-026-01439-z",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
19
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41531-026-01380-1 [code]
Identifying maximal beta power from directional subthalamic local field potentials in Parkinson's disease.
Journal: NPJ Parkinson's disease
In common: MNE-Python, pandas, SciPy, 2 other tools, Parkinson's, clinical / translational, 2 authors
[2] doi:10.1038/s41591-026-04434-2 [code]
Adaptive deep brain stimulation for dynamic gait control in Parkinson's disease: a randomized feasibility trial.
Journal: Nature medicine
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, pandas, 3 other tools, Parkinson's, clinical / translational, 4 references
[3] doi:10.1038/s41591-026-04432-4 [code]
Activity-dependent adaptive deep brain stimulation improves gait in Parkinson's disease.
Journal: Nature medicine
In common: MNE-Python, SciPy, Matplotlib, 1 other tool, Parkinson's, clinical / translational, 4 references
[4] doi:10.1002/ana.78206 [code]
Multimodal Image Guidance in Subthalamic Deep Brain Stimulation for Parkinson's Disease.
Journal: Annals of neurology
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, pandas, 3 other tools, Parkinson's, clinical / translational, author Andrea A Kühn
[5] doi:10.1111/ene.70678 [code]
Who Falls After a Stroke? Evidence From a Prospective Stroke Cohort.
Journal: European journal of neurology
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, pandas, 3 other tools, clinical / translational, author Andrea A Kühn
[6] doi:10.1016/j.celrep.2026.117404 [code]
Action and rest tremor map to distinct networks within the primary motor cortex.
Journal: Cell reports
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, pandas, 3 other tools, author Andrea A Kühn
[7] doi:10.1126/sciadv.aea3919 [code]
Hierarchical brain dynamics supporting visual perceptual transitions.
Journal: Science advances
In common: MNE-Python, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 4 other tools, 1 reference
[8] doi:10.64898/2026.03.06.710026 [code]
Distinct beta burst motifs exhibit opposing error relationships during motor adaptation
Journal: bioRxiv (preprint)
In common: MNE-Python, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 4 other tools, 1 reference
[9] doi:10.1038/s41597-025-06397-4 [code]
MEG-SCANS - A comprehensive magnetoencephalography speech dataset with Stories, Chirps and Noisy Sentences
Journal: n/a
In common: MNE-Python, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 4 other tools, 1 reference
[10] doi:10.1162/imag.a.1040 [code]
Novel 4 He-OPMs support waveform-specific beta burst analysis comparable to SQUID-MEG
Journal: n/a
In common: MNE-Python, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 4 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.