<jupyter_start><jupyter_text>![SegmentLocal](../../assets/images/Logo2.png)  Simulating Data*Neural Time Series Data*  PrerequisitesFor this chapter, you should be familiar with the following concepts and techniques:* Basic Python programming* Basic Math. **(recap your skills in Linea Algebra, Sine Waves and Euler's Formula)**  Scope of this tutorialIn this tutorial, you will learn the conceptual, mathematical, and implementational (via python programming) basis of time- and time-frequency-analysis of EEG recordings. Alternating between theoretical background knowledge and practical exercises you will learn the basics of how EEG Data is recorded, preprocessed and analysed. We will however only cover the **fundamental basics of EEG analysis**; but with this, you will then be prepared to dig deeper into the endless opportunities of Neural Time Series analysis.General remark: In order to make the most out of the exercises, we highly recommend you to exploit the interactive character of these notebooks; play around with the parameters and see how this affects your results... (e.g. what happens if you use a lower transition bandwidth in your filter or if you don't extract the abs() values of the Fourier Transform?) This way, it is much easier to really understand what you are doing!  1. The Math behind EEG signalsBefore we are ready to work with real EEG data, we will first create artificial signals. This makes it much easier to understand the *maths* behind EEG signals, which in return will help you to understand the following  analysis steps a lot better.Step by step, we will make our signal more complex until it approximates *'real'* EEG data. In the next section of this chapter, we will then start to use this knowledge in order to analyse EEG Data, recorded by a Neurobiopsychology research group of our institute.For the following exercises, we will use the signal processing toolbox from scipy. [This link](https://docs.scipy.org/doc/scipy/reference/signal.html) leads you to the documentation of the toolbox, where you can find almost all the functions that you need to solve the following tasks. Whenever you are working with a toolbox, I highly recommend to take some time to explore the corresponding documentation. It helps you to make the most of all the tools it supplies!<jupyter_code>!pip install mne 
#run this if mne is not installed on the colab kernel
import matplotlib.pyplot as plt
import numpy as np
import random
from scipy import signal 
#dont name any of your variables signal otherwise the package won't work!<jupyter_output><empty_output><jupyter_text>Simple signals One of the key concepts you need in order to understand the maths behind oscillatory signals (like neural signals in EEG) is the **sine wave**. Here you can find a short overview of the parameters that define a sine wave (more details have been covered in the video "Analysing Neural Time Series Data / EEG Intro" in Chapter 6.1. on studIP).  ![SineUrl](https://media.giphy.com/media/U6prF59vkch44/giphy.gif "Sine") ![SegmentLocal2](../../assets/images/sinewave.png) With the parameters ```amplitude```,```frequency``` and ```phase```($\theta$), a sine wave can be described with the following formula:$$Asin(2*\pi ft + \theta)$$  1.1 Simple SinewaveWith this information, we are now ready to create a simple signal as a combination of two sinusoids. For this:  - Define a time scale of 1 second, i.e. 1000ms- Create THREE sinewaves with a length of 1sec: one with a frequency of 10Hz, 15Hz, and 20Hz (*for simplicity, we will for now ignore amplitude and phase, they will be used in the next step though*)- Add them together to create your first simple signal- Create a plot for each of the sinusoids<jupyter_code>t = np.linspace(0, 1, 1000, False)  # 1 second
...

# plot all three figures
fig, ax = plt.subplots(4, 1, sharex=True, sharey=True, figsize=(16,8))

# Plot each graph
...

plt.show()<jupyter_output><empty_output><jupyter_text>1.2 More complex signalAs a next step, we want to achieve something more 'complex'.For this, we select a list of frequencies, that we want our signal to be composed of, and define their amplitudes and phases. The exact values that you should use for this are already predefined; but play around with them and see how your results change!With the help of these parameters:- Create a new, a bit more complex signal by combining the resulting sinusoids (you should get 6 sinusoids with the respective ```freq```, ```amplit``` and ```phase```)- To make it more realistic, create some random Gaussian noise with the same length and add it to your signal- Then plot both, the clean and the noisy signalThe Nyquist Sampling Theorem states: In order to prevent distortions of the underyling information, the minimum sampling frequency of a signal should be double the frequency of its highest frequency component (which is in our case 60Hz).We will define it a little bit higher for our filter to properly work. But you can of course change it and see how this affects your plots (especially your filter in exercise 2.2).<jupyter_code># set parameters:
srate = 1000 # define sampling rate 
nyq = srate/2 #nyquist frequency
freq = [3, 10, 5, 15, 35, 60] # define a list of frequencies
amplit = [5, 15, 10, 5, 7, 1] # define their amplitudes
phase = [np.pi/7, np.pi/8, np.pi, np.pi/2, -np.pi/4, np.pi/3] # and their respective phases
# 1. create signal
t = np.linspace(0, 1, 1000, False)  # 1 second
sig = []

# HINT: use a for loop
...

# 2. add some random noise
# info: the third paramter defines the size of your noise HINT: use np.random.normal
noise = ...

signal_final = ...
signal_noisy = ...

# 3. plot both figures (signal with and without noise)
fig, axs = plt.subplots(2, 1, sharex=True, sharey=True, figsize=(15,8))

# add big shared axes, hide frame to share ylabel between subplots
fig.add_subplot(111, frameon=False)
plt.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
plt.xlabel('Time [s]',fontsize = 12)
plt.ylabel('Amplitude',fontsize = 12)

# Plot each graph
...
plt.show()<jupyter_output><empty_output><jupyter_text>In reality, your EEG signal is roughly based on the same components: it typically contains a mixture of simultaneous neural oscillations at different frequencies plus some noise. This noise can be non-neural (caused by line noise or muscle activity); but also neural oscillations that are not of your interest can be considered as 'noise'. In order to be able to do your analysis as "clean" as possible, you want to isolate only the part of the signal that you are interested in, thereby increasing the **signal-to-noise-ratio (SNR)** of your signal. A way to do this, is by **filtering** your data - this will be the focus of the following exercises.  2. How to get rid of unwanted noise?  The Fourier Transform Before telling you more about how EEG data can be filtered, you need to first learn about the **Fourier Transform (FT)**, which is a really useful and important mathematical tool to analyse EEG time series data (and any kind of time series data in general); with its help we can separate the different frequency components that compose our signal and thus get rid of unwanted frequencies. To get a more intuitive explanation of Fourier Transform watch this video by [3Blue1Brown](https://www.youtube.com/watch?v=spUNpyF58BY&ab_channel=3Blue1Brown).  → If you want to have a more thorough and more EEG-based intro, [check out this vdeo](https://www.youtube.com/watch?v=d1Yj_7ti_IU&list=PLn0OLiymPak3lrIErlYVnIc3pGTwgt_ml&index=3) by <NAME>. Or [watch the whole playlist](https://www.youtube.com/playlist?list=PLn0OLiymPak3lrIErlYVnIc3pGTwgt_ml) for even more detailed explanations.  2.1 Extracting the frequency spectrum with the FFT Now we are ready to apply the **fast Fourier Transform** [```fft.fft()```](https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html) on our signals, in order to get its frequency spectrum. Since we created the signal on our own, we can check whether it shows all the frequencies that we used to compose it. * We will have to compute the frequency spectrum several times, therefore it is useful to write a proper function that can compute the fourier transform of any signal and already selects the absolute part of the Fourier Transform. The explanation for why we only use the absolute values of the FT gets really mathematical; but if you are curious, I can highly recommend [this video!](https://www.youtube.com/watch?v=Nupda1rm01Y).* Then apply and plot the FFT of the noisy signal from 1.2. You will see a lot more, if you limit your x-axis to 100, since we are not interested in higher frequencies anyway.If you look at the amplitudes, you will realize that they are half the amplitudes that we predefined when creating the signal. This happens because we are only taking the absolute values of the FT-frequencies and the amplitudes for the negative frequencies have been "removed".<jupyter_code>def getFT(sig):

    # compute fft
    FFT = np.fft.fft(sig)
    FFT = np.abs(FFT)/len(sig) #normalise the frequencies to the number of time-points
    
    return FFT
# compute and plot FFT of the noisy signal

frequencies = getFT(...)
N = int(len(frequencies)/2)

fig, ax = plt.subplots(figsize=(8,8))
ax.plot(frequencies[:N]) #What happens when you plot all frequencies? Try it!!!
plt.suptitle('Fourier Transform of the Signal')
ax.set(
    title='Frequencies {} plus random Gaussian noise'.format(freq),
    xlim=(0,100),
    xlabel='Frequency (Hz)',
    ylabel='Amplitude'
)

plt.show()
print(N)<jupyter_output><empty_output><jupyter_text>Filtering EEG Data Now that we have recovered the components of our signal with Fourier Transform, we will have a short look into how EEG data is filtered,in order to remove the noise. This knowledge will be important for the second half of the Notebook.  2.2 Filtering in the time-domain vs. filtering in the frequency-domain In this part, we will see with some practical examples how we can filter our data in two different ways. Usually, you only filter in the frequency domain since this is computationally a lot faster. Yet, it is really useful to learn about both procedures in order to better understand the concept of filtering in general. In the video above you already got a first impression of filtering in the frequency domain. In order to better understand its time-domain equivalent, you need to first learn about the process of convolution, i.e. the (mathematical) procedure of applying your filter to your data in the time domain. This should however not be entirely new to you, you will find some similarities to the procedure behind the Fourier Transform:Convolution: Convolution is used to isolate frequency-band-specific activity and to localize that frequency-band-specific activity in time.  This is done by convolving wavelets— i.e. time-limited sine waves—with EEG data. As the wavelet (i.e. the convolution kernel) is dragged along the EEG data (the convolution signal): it reveals when and to what extent the EEG data contain features that look like the wavelet. When convolution is repeated on the same EEG data using wavelets of different frequencies, a time-frequency representation can be formed."   (<NAME>, "Analyzing Neural Time Series Data: Theory and Practice"→ If you want a more thorough and more visual explanation of convolution, I can highly recommend [this video](https://www.youtube.com/watch?v=9Hk-RAIzOaw) by <NAME>.Convolution theorem: Convolution in the time domain is the same as multiplication in the frequency domain.  ![SegmentLocal](../../assets/images/Convolution_Theorem.png)  2.2.1 Filter in the time domain According to the figure above, in order to filter our signal in the time domain, we need a windowed sinewave as a filter-kernel. The windowing helps to obtain **temporally localized frequency information**. We then convolve this wavelet with our signal, extracting the frequency bands that we want to work with. - First define your pass-band as 25Hz. Ideally everything above this frequency is filtered out; in reality however, we need a transition band of about 10 Hz, or a region between the pass-frequency ```f_p``` and stop-frequency ```f_s```. In this range, frequencies are only attenuated instead completely excluded. This is necessary in order to account for the trade-off between precision in the frequency-domain and precision in the time-domain. - Next, we define the gains at each frequency band: everything outside 0 and our pass-band of 25Hz should be attenuated, i.e. have a gain close to 0.- Using the function [```firwin2()```](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.firwin2.html) of the signal package and the parameters from above, we can now construct our filter-kernel ```h_win``` (the result should be a wavelet with a length/duration of 0.6 seconds) - Plot your kernel as well as its frequency spectrum. It should look like a step-function, that assigns a gain of 1 to all frequencies in our pass-band between 0 - 25Hz.Tip: Play around with your the parameters of your filter (e.g. the filter's duration, its transition bandwidth or its stop- and passband, the sampling rate etc.) and see how the plots change. You can then also proceed with the whole filtering process and check out what different filters (with different paramters) do to your data. This way, you can properly understand how the different parameters are finally affecting you data.<jupyter_code># Create a Low-pass Filter: Windowed 10-Hz transition

# 1. Define Filtering parameteres
filter_duration = 0.8

n = int(round(srate * filter_duration)+1) #odd number of filter coefficients for linear phase
f_p = 25. # define passband
trans_bandwidth = 10.
f_s = f_p + trans_bandwidth  # stopband = 35 Hz
print(f_s)
# define gains of each frequency band
freq = [0., f_p, f_s, nyq]
gain = [1., 1., 0., 0.]

# 2. Compute filter graph
h_win = signal.firwin2(..., nyq=nyq)
# 3. Compute freqeuncy spectrum of the filter
frequencies = getFT(...) 

# 4. Plot filter in time and in frequency domain
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,8))
plt.subplots_adjust(hspace=0.5)
time =  np.linspace(-1, 1, len(h_win))

...

plt.show()

# use inbuilt mne-function to plot filter characteristics (a lot more detailed)
import mne
flim = (1., nyq) # limits for plotting
mne.viz.plot_filter(h_win, srate, freq, gain, 'Windowed {} Hz transition ({}s filter-duration)'.format(trans_bandwidth, filter_duration),flim=flim, compensate=True)<jupyter_output><empty_output><jupyter_text>Now we are ready to convolve our signal with our self-constructed FIR filter ```h_win```. - For this, we use the [```convolve()```](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve.html) function of the *signal* package. - Then plot both, the filtered and the unfiltered signal in order to see the effect of your filter.- Aftwards, we want to test in the frequency-spectrum of our signal whether our filter successfully attenuated the frequency-components above 25Hz. For this: compute and plot the FT of both, the filtered and the unfiltered signal. - In order to compare which filtering procedure is faster, record the computation-time of the time-domain convolution with the help of the magic function [```%timeit```](https://docs.python.org/2/library/timeit.html) (you can write an extra line for this, where you again perform your convolution).<jupyter_code># 1. Convolve signal with the filter
conv_time = signal.convolve(..., mode='same')

# and Calculate computation time
%timeit signal.convolve(..., mode='same')
# 2. Plot filtered and unfiltered signal
...
# 3. Compute and plot frequency spectrum of the filtered and the unfiltered signal 
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,8))
plt.subplots_adjust(hspace=0.5)

...

plt.show()<jupyter_output><empty_output><jupyter_text>2.2.2 Filter in the frequency domainFiltering in the frequency domain is computationally much faster and easier. According to the convolution theorem (see above):- Multiply the frequency-spectrum of our filter-kernel with the frequency-spectrum of our signal.- In order to compare the filtered and the unfiltered signal, first compute the inverse Fourier-Transform of your filtering result with [```fft.ifft```](https://numpy.org/doc/stable/reference/generated/numpy.fft.ifft.html) (in order to translate it to the time domain) and then plot both signals, the unfiltered and filtered one, in one plot.  Note: So far, every time we applied the fourier transform (FT) to our signal, we only used the absolute values of the FT-result, because this was what we were interested in. To visualize what that means, just plot the FT of any of our signals with and without the abs()-function. For the inverse FT to work porperly however, we need the "whole" result of the FT, which is why we omit the abs() function this time.- In a second plot, compare your result from filtering in the frequency domain with your convolution result from before in the time domain (from 2.2 a). According to relationship between frequency and time domain, both curves should look exactly the same!- In order to compare which filtering procedure is faster, again record the computation-time of the frequency-domain filtering with the help of the magic function [```%timeit```](https://docs.python.org/2/library/timeit.html). Compare the result to the computation time of the time-domain convolution. Which one is faster?<jupyter_code># 1. Compute lengths of the result
# in order to make the inverse FFT return the correct number of time points:
# we need to make sure to compute the FFTs of the signal and the kernel using 
# the appropriate number of time points. In other words:
# the length of the signal (=srate = 1000) plus the length of the kernel (= 801)  minus one (= 1800)
# Afterwards the result has to be trimmed to its orignal length again (step 5)
n_signal = ...
n_kernel = ...
print(n_signal, n_kernel)
nconv = n_signal + n_kernel -1
halfk = np.floor(n_kernel/2)

# 2. Compute FT of the kernel and the signal
h_winX = np.fft.fft(...)
signalX = np.fft.fft(...)

# 3. Multiply frequecies
result_frequencydomain = h_winX*signalX

# 4. Compute inverse FT (convert frequency-domain to the time-domain)
result_timedomain = np.fft.ifft(result_frequencydomain)
print(len(result_timedomain))

# 5. Cut the signal to original length
result_timedomain = result_timedomain[int(halfk):-int(halfk)]

# 6. Plot both signals (unfiltered and filtered) in one plot
...
plt.legend(bbox_to_anchor=[1.2, 1.3], ncol=2)

# 7. Plot results of filtering in the frequency domain and filtering in the time domain

fig, ax = plt.subplots(3,1, figsize=(15,8), sharey=True)

...


plt.subplots_adjust(hspace=0.2)
fig.add_subplot(111, frameon=False)
plt.suptitle('Compare results of filtering in the frequency domain and filtering in the time domain') 
plt.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
plt.xlabel('Time [ms]',fontsize = 12)
plt.ylabel('Amplitude',fontsize = 12)
plt.legend();
# 8. calculate computation time
filtering_fd_time = get_ipython().run_line_magic('timeit', "conv_fd = h_winX*signalX")<jupyter_output>3.06 µs ± 46 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)<jupyter_text>Bonus Exercise: Time-Frequency Analysis FT alone does not describe the signal perfectly. For non-stationary signals (like EEG), we are interested in the evoked response of the brain. The result of the simple FT alone will not show us that.  Hence, we rely on Time-Frequency Analysis in order to understand the temporal structure of the different frequencies in the signal. **Spectrograms** will do the trick! After applying a specific (time-resolved) version of the FT, they shows us how much of each frequency compoment was present at a specific time point. [This section of one of Cohen's videos about the Fast Fourier Transform](https://youtu.be/T9x2rvdhaIE?t=118), nicely explains and visualizes how a spectogram is computed.- Plot the power spectogram of the noisy signal using the function [```plt.specgram()```](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.specgram.html).- Also compare the results with the [stft](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.stft.html) method<jupyter_code># Plot simple spectogram of the noisy signal
...<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07841504628495861, 'fraction_numerical': 0.013126316305039918, 'mean_word_length': 4.701759285115889, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 35, 'https://': 16, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '22828', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 5894, 'n_tokens_neox': 5462, 'n_words': 2811, 'path': 'Tutorials/Day-04/Data_simulation.ipynb', 'repo_name': 'ooezenoglu/ACC_SS2021', 'size': '20436'}