Coverage for brainbox/lfp.py: 0%

1# -*- coding: utf-8 -*- 

2""" 

3Functions to analyse LFP signals. 

4 

5@author: Guido Meijer 

6Created on Fri Mar 13 14:57:53 2020 

7""" 

8 

9from scipy.signal import welch, csd, filtfilt, butter 

10import numpy as np 

11 

12 

13def butter_filter(signal, highpass_freq=None, lowpass_freq=None, order=4, fs=2500): 

14 

15 # The filter type is determined according to the values of cut-off frequencies 

16 Fn = fs / 2. 

17 if lowpass_freq and highpass_freq: 

18 if highpass_freq < lowpass_freq: 

19 Wn = (highpass_freq / Fn, lowpass_freq / Fn) 

20 btype = 'bandpass' 

21 else: 

22 Wn = (lowpass_freq / Fn, highpass_freq / Fn) 

23 btype = 'bandstop' 

24 elif lowpass_freq: 

25 Wn = lowpass_freq / Fn 

26 btype = 'lowpass' 

27 elif highpass_freq: 

28 Wn = highpass_freq / Fn 

29 btype = 'highpass' 

30 else: 

31 raise ValueError("Either highpass_freq or lowpass_freq must be given") 

32 

33 # Filter signal 

34 b, a = butter(order, Wn, btype=btype, output='ba') 

35 filtered_data = filtfilt(b=b, a=a, x=signal, axis=1) 

36 

37 return filtered_data 

38 

39 

40def power_spectrum(signal, fs=2500, segment_length=0.5, segment_overlap=0.5, scaling='density'): 

41 """ 

42 Calculate the power spectrum of an LFP signal 

43 

44 Parameters 

45 ---------- 

46 signal : 2D array 

47 LFP signal from different channels in V with dimensions (channels X samples) 

48 fs : int 

49 Sampling frequency 

50 segment_length : float 

51 Length of the segments for which the spectral density is calcualted in seconds 

52 segment_overlap : float 

53 Fraction of overlap between the segments represented as a float number between 0 (no 

54 overlap) and 1 (complete overlap) 

55 

56 Returns 

57 ---------- 

58 freqs : 1D array 

59 Frequencies for which the spectral density is calculated 

60 psd : 2D array 

61 Power spectrum in V^2 with dimensions (channels X frequencies) 

62 

63 """ 

64 

65 # Transform segment from seconds to samples 

66 segment_samples = int(fs * segment_length) 

67 overlap_samples = int(segment_overlap * segment_samples) 

68 

69 # Calculate power spectrum 

70 freqs, psd = welch(signal, fs=fs, nperseg=segment_samples, noverlap=overlap_samples, 

71 scaling=scaling) 

72 return freqs, psd 

73 

74 

75def coherence(signal_a, signal_b, fs=2500, segment_length=1, segment_overlap=0.5): 

76 """ 

77 Calculate the coherence between two LFP signals 

78 

79 Parameters 

80 ---------- 

81 signal_a : 1D array 

82 LFP signal from different channels with dimensions (channels X samples) 

83 fs : int 

84 Sampling frequency 

85 segment_length : float 

86 Length of the segments for which the spectral density is calcualted in seconds 

87 segment_overlap : float 

88 Fraction of overlap between the segments represented as a float number between 0 (no 

89 overlap) and 1 (complete overlap) 

90 

91 Returns 

92 ---------- 

93 freqs : 1D array 

94 Frequencies for which the coherence is calculated 

95 coherence : 1D array 

96 Coherence takes a value between 0 and 1, with 0 or 1 representing no or perfect coherence, 

97 respectively 

98 phase_lag : 1D array 

99 Estimate of phase lag in radian between the input time series for each frequency 

100 

101 """ 

102 

103 # Transform segment from seconds to samples 

104 segment_samples = int(fs * segment_length) 

105 overlap_samples = int(segment_overlap * segment_samples) 

106 

107 # Calculate coherence 

108 freqs, Pxx = welch(signal_a, fs=fs, nperseg=segment_samples, noverlap=overlap_samples) 

109 _, Pyy = welch(signal_b, fs=fs, nperseg=segment_samples, noverlap=overlap_samples) 

110 _, Pxy = csd(signal_a, signal_b, fs=fs, nperseg=segment_samples, noverlap=overlap_samples) 

111 coherence = np.abs(Pxy) ** 2 / (Pxx * Pyy) 

112 phase_lag = np.angle(Pxy) 

113 

114 return freqs, coherence, phase_lag