<jupyter_start><jupyter_text>T11 Example project with tips and tricks  Hint and tips are provided with yellow background.  Project : Extract local field potential (LFP) from extracellular recording in the hippocampus You are provided with a raw extracellular recording from the hippocampus of a awake behaving rat. The trace has been recorded with a silicon probe (32 channel) and you will be analyzing the signal from a single channel (the probe delivers in total 32 of such recordings). The signal contains the local field potential (LFP) as well as spiking signals. Since both occur at different frequencies, it is possible to separate them through appropriate signal filtering using spectral analysis. Your task is to extract frequency content of the LFP from the raw trace. The frequency content of the LFP trace can then be further analyzed and be linked to specific behaviors.  Objective of the projectStudy the frequencies contained in the LFP signal. For that purpose, calculate the spectrogram showing the frequency content as a function of time. The frequency range of interest ranges from 0 to 100 Hz. Characterize which frequencies are contained in the recording, and whether these signals change over time. Which frequencies show prominent dynamics during the recording period, i.e., which frequencies show dynamics changes in spectral power over time?   Using the information provided, try to learn as much as possible about the scientific and technical background of the project. * have a in depth look at the publication provided with the project (if any), here : [<NAME> et al. (2017) Nature Neurosciene](https://www.nature.com/articles/nn.4637)* search for project keywords, here : "hippocampus", "silcon probe", "local field potential", "spectral analysis", "spectrogram"* look back at the lecture and turorial which dealt with the technique to be used in the project, here : "Course 6 : Spectral analysis of analog signals"   Details about the recording The data has been recorded in the lab of [<NAME>](https://girardeaulab.org/). The entire trace is 6.3 hours long with running on track at the beginning, middle and end and homecage inbetween with sleep and quiet wakefulness. The electrode has been placed in the hippocampus in the CA1 pyramidal cell layer. The sampling rate of the recording is 20 kHz (or 20 000 Hz). To avoid issues with handling the large data file, a part of the entire recording was extracted for this project.   The data Find below the code on how to load the data.  Try to first understand in detail the structure of the data. What is contained in the data structure and what are the properties of the data.<jupyter_code>import numpy as np    
import matplotlib.pyplot as plt
# if you use qt (instead of inline what we used in the course) you can zoom in on the data which is useful for large traces
%matplotlib qt 

data = np.load('data-for-project_3.npy')

print(np.shape(data))
print(data[0])
print(data[1])

# look in more detail at the time array
samplingRate = 1./np.mean(np.diff(data[0]))

print('The sampling rate is : ', samplingRate, ' Hz')

print('The total duration of the recording is ', data[0,-1], ' s. Or ', data[0,-1]/60 ,' min')

# look in more detail at the data array

print('max :',np.max(data[1]))
print('min :',np.min(data[1]))
print('mean : ', np.mean(data[1]))<jupyter_output>(2, 36000000)
[0.00000000e+00 5.00000006e-05 1.00000000e-04 ... 1.79999985e+03
 1.79999990e+03 1.79999995e+03]
[898. 882. 882. ... 502. 518. 522.]
The sampling rate is :  19999.999999999996  Hz
The total duration of the recording is  1799.9999500000004  s. Or  29.999999166666672  min
max : 19146.0
min : -16198.0
mean :  24.52305088888889<jupyter_text>Plotting Plot the raw data!As always, plotting the data provide immediate and useful insights. Generally, plot data and results helps to develop a better understanding. I would recommend plotting in an external window, since this allows to zoom into the data. This is in particular useful for large data-sets. Don't forget to include figures in the jupyter-notebook for the project presentation and for illustrating the results obtained.  * zoom into the data and explore features at many different time points, here : check out the interval [70,90] sec * try to get a first rough estimate about the project question through visual inspection, here : what is the frequency present in the signal  * if data-set is too large, plot and inspect smaller sections of the data  * figures can be implemented with the code : ![Data Figure](Figure_data.png)<jupyter_code>plt.plot(data[0],data[1]) # plot the first 1000 entries of the time series
plt.xlabel('time (s)')
plt.show()<jupyter_output><empty_output><jupyter_text>AnalysisUsing the course material, try to understand the analysis to perform and how to implement it. * Pick the right lecture and read up about the analyis method. Maybe read more extensive explanations on the internet about the analysis, here : [Spectral analysis on Wikipedia, for exmaple](https://en.wikipedia.org/wiki/Spectral_density_estimation)* Use example code snippets from the tutorial to start programming, here : "T06 : Calculating the spectrogram of a signal"* Go through the code and try to understand the code. * You could also start from a code snippet you find on the internet. There are many possible ways to realize similar analysis. The difficulty is to find what suits best your problem.<jupyter_code>import scipy.signal as sg


# constructing the signal #################
fs = 1000                          #1
time = np.linspace(0,2,fs*2)       #2
y=sg.chirp(time,400,1,0,'linear')  #3


# computing the spetrogram ################
windLength = 128              #4
overl = windLength-1          #5
wind=np.hanning(windLength)   #6
#wind=np.kaiser(windLength,0) #7
f, tt, Sxx =sg.spectrogram(y,fs,wind,len(wind),overl) #8

# producing the figure ################
fig=plt.figure(figsize=(10,10))  #9

ax0 = plt.subplot(2,1,1)          #10
ax0.plot(time,y)                  #11
ax0.set_xlabel('time (s)')        #12
ax0.set_ylabel('signal')          #13

ax1 = plt.subplot(2,1,2)          #14


im = ax1.pcolormesh(tt,f,Sxx)     #15
fig.colorbar(im, ax=ax1)          #16  
ax1.set_xlabel('time (s)')         #17
ax1.set_ylabel('freqeuncy (Hz)')            #18

plt.show()                        #19<jupyter_output><empty_output><jupyter_text>Adapt the code for the project. In other words, start from the tutorial code and adapt it  * Change the values and fill in the project specifics. * It is useful and good practise to separate **computations**: the actual implementation of the analysis and **visualization** : the figure generation.  **Computation**<jupyter_code>import scipy.signal as sg


# constructing the signal #################
fs = 20000                         

# computing the spetrogram ################
windLength = 128              #4
overl = windLength-1          #5
wind=np.hanning(windLength)   #6
f, tt, Sxx =sg.spectrogram(data[1],fs,wind,len(wind),overl) #8
#f, tt, Sxx =sg.spectrogram(data[1], fs, window='hanning',nperseg=1024, noverlap=1024 - 100,detrend=False, scaling='spectrum')<jupyter_output><empty_output><jupyter_text>**Figure generation**<jupyter_code># producing the figure ################
fig=plt.figure(figsize=(10,10))  #9

ax1 = plt.subplot(1,1,1)          #14

im = ax1.pcolormesh(tt,f,Sxx)     #15
fig.colorbar(im, ax=ax1)          #16  
ax1.set_xlabel('time (s)')         #17
ax1.set_ylabel('freqeuncy (Hz)')            #18

plt.show()                        #19<jupyter_output><empty_output><jupyter_text>The above code crashes the jupyter-notebook! It also aborts when executing the code in `ipython` for example. There exits probably a memory problem : the execution of the code requires more RAM than available on the computer. **Important lesson :** Programming is trial and error. Most first tries don't always provide the desired outcome. Most imporantly, don't panic. This is normal!Your tutor might have warned you about possible pitfalls. If not, play around with the code for a while and if no solution is found, reach out the the tutor.  Alternative method to analyze the data We could split the signal is segments and calculate the frequency content sequentially. We start by chopping up the signal into slices of 1024 samples, each slice overlapping the previous by 100 samples. The resulting slices array contains one slice per row. We found a `scikit-learn` function online to perform that : [Frequency and the Fast Fourier Transform](https://www.oreilly.com/library/view/elegant-scipy/9781491922927/ch04.html)<jupyter_code>from skimage import util

windowLength = 1024

slices = util.view_as_windows(data[1], window_shape=(windowLength,), step=100)

print(np.shape(data[1]))
print(np.shape(slices))
print(windowLength)

nSlices = len(slices)
print(nSlices)<jupyter_output>(36000000,)
(359990, 1024)
1024
359990<jupyter_text>We have another code snipplet from tutorial 6 where we calculate the spetral content using FFT in numpy. In most cases, it is more efficient and instructive to perform the analysis first on a subset of the data. * Optimize the code and visualization until the desired output is generated. Once the output is plausible and correct, move on to apply the analysis to the entire data.<jupyter_code>hanningWindow = np.hanning(windowLength)
fs = 20000
exampleSlice = slices[1]*hanningWindow

nyquist = fs/2.                                  #10
fSpaceSignal = np.fft.fft(exampleSlice)               #11
fBase = np.linspace(0,nyquist,np.floor(len(exampleSlice)/2)+1) #12
halfTheSignal = fSpaceSignal[:len(fBase)]        #13
complexConjugate = np.conj(halfTheSignal)        #14
powe = np.abs(halfTheSignal*complexConjugate)            #15

plt.plot(fBase,powe/max(powe))
plt.show()
print(fBase[:20])<jupyter_output>[  0.       19.53125  39.0625   58.59375  78.125    97.65625 117.1875
 136.71875 156.25    175.78125 195.3125  214.84375 234.375   253.90625
 273.4375  292.96875 312.5     332.03125 351.5625  371.09375]<jupyter_text>Here in particular, we need to change the resolution of the frequency axis. This can be done by changing the window length.<jupyter_code>windowLength = 2**16

slices = util.view_as_windows(data[1], window_shape=(windowLength,), step=100)

hanningWindow = np.hanning(windowLength)

fs = 20000
exampleSlice = slices[1]*hanningWindow
nyquist = fs/2.                                  #10
fSpaceSignal = np.fft.fft(exampleSlice)               #11
fBase = np.linspace(0,nyquist,np.floor(len(exampleSlice)/2)+1) #12
halfTheSignal = fSpaceSignal[:len(fBase)]        #13
complexConjugate = np.conj(halfTheSignal)        #14
powe = np.abs(halfTheSignal*complexConjugate)            #15

print(fBase)
print(powe)<jupyter_output>[0.00000000e+00 3.05175781e-01 6.10351562e-01 ... 9.99938965e+03
 9.99969482e+03 1.00000000e+04]
[1.03231620e+13 3.49448185e+13 3.52741284e+13 ... 2.70800143e+06
 3.37002531e+06 5.73780545e+06]<jupyter_text>Further optimize We are interested in the frequency range up to 100 Hz, so lets restrict the output to that.<jupyter_code>mask = fBase < 100
plt.plot(fBase[mask],powe[mask]/max(powe[mask]))
plt.show()<jupyter_output><empty_output><jupyter_text>The specific unit used for the ratio is the decibel, `20*log10` (amplitude ratio).<jupyter_code>S = 20 * np.log10(powe[mask] / np.max(powe[mask]))

plt.plot(fBase[mask],S)
plt.show()<jupyter_output><empty_output><jupyter_text>Extend the analysis to the entire data-setUse the above developped code and apply it to the entire data-set. Sometimes, the computation can take minutes or even hours, when the data-set is large or the computations heavy. Be patient. You can try to let the analysis run over night. * Copy all relevant code in one window. That way you avoid surprises of a stray variable which has been initialized with a strange value somewhere else.  For the above problem, let's write a loop which calculates spetral power for all slices.<jupyter_code>import pdb
fs = 20000
windowLength = 2**16
import scipy.signal as sg 

slices = util.view_as_windows(data[1], window_shape=(windowLength,), step=int(windowLength*0.9))

nSlices = len(slices)

hanningWindow = np.hanning(windowLength)
nyquist = fs/2.  
mask = fBase < 100

spectra = []

for i in range(nSlices):
    exampleSlice = slices[i]*hanningWindow                               #10
    fSpaceSignal = np.fft.fft(exampleSlice)               #11
    fBase = np.linspace(0,nyquist,np.floor(len(exampleSlice)/2)+1) #12
    halfTheSignal = fSpaceSignal[:len(fBase)]        #13
    complexConjugate = np.conj(halfTheSignal)        #14
    powe = np.abs(halfTheSignal*complexConjugate)            #15
    #S = 20 * np.log10(powe[mask] / np.max(powe[mask]))
    #S = powe[mask] / np.max(powe[mask])
    spectra.append(powe[mask])

frequencies = fBase[mask]
spectrum = np.asarray(spectra)
psd = spectrum/np.max(spectrum)
np.max(spectrum)<jupyter_output><empty_output><jupyter_text>Figure generation Show your results in figure with captions and labels.  Producing good representations of results in form of figures takes time and effort. Try to find the visualization which best conveys your result. Certain measures have standard representations (e.g. frequency content is shown in spectrograms, spikes in raster plots, etc.).<jupyter_code>plt.imshow(psd.T, origin='lower', cmap='viridis',
          extent=(0, 30*60, 0, 100),aspect='auto',interpolation='nearest')
plt.colorbar()
plt.xlabel('time (sec)')
plt.ylabel('frequency (Hz)')<jupyter_output><empty_output><jupyter_text>Refine and interpret results The specific objective for this project was to find frequency ranges which exhibit dynamics over time. Here, it looks like as if theta power is present only for certain time intervals. Let's extract the mean power in the theta range and display it over time. Going an additional step after obtaining the main result provides often great insight, doesn't take much effort and increases the quality/impact of your work.<jupyter_code>plt.imshow(psd.T, origin='lower', cmap='viridis',
          extent=(0, 30*60, 0, 100),aspect='auto',interpolation='nearest')
plt.axhline(y=5,ls='--',color='black')
plt.axhline(y=10,ls='--',color='black')
plt.colorbar()
plt.xlabel('time (sec)')
plt.ylabel('frequency (Hz)')
thetaMask = (frequencies>6)&(frequencies<9)

print(frequencies[thetaMask])

timeEntries = np.shape(psd.T)[1]
time = np.linspace(0,30,timeEntries)

plt.plot(time,np.mean((psd.T)[thetaMask],axis=0))
plt.xlabel('time (min)')
plt.ylabel('mean theta power (dB)')<jupyter_output>[6.10351562 6.40869141 6.71386719 7.01904297 7.32421875 7.62939453
 7.93457031 8.23974609 8.54492188 8.85009766]
{'fraction_non_alphanumeric': 0.09965916837082481, 'fraction_numerical': 0.06039536468984322, 'mean_word_length': 4.325226860254084, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 51, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 8, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '54581', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 5106, 'n_tokens_neox': 4415, 'n_words': 1646, 'path': 'tutorials/T11_Project-LFP_from_raw_ephys_trace.ipynb', 'repo_name': 'mgraupe/DataSciPy2021', 'size': '14696'}