Power Spectral Desities of EEG signals

This algorithm allows you to calculate the power spectral density of an EEG signal in a very simple way.

The present algorithm allows graphing the power spectra of the alpha and beta brain rhythms, which were acquired during a Motor Imagery (MI) experiment. This algorithm is structured in three phases that are; pre-processing, conformation of the alpha and beta oscillations, and estimation of the power spectral density for said rhythms.

During the acquisition of the EEG signals, the MI stages were constituted by ten trials, each one of them is in turn made up of twelve runs, thus, a run is structured by three phases that are: Relax, Quietude and MI,each one of three seconds of duration.

In the algorithm the traces acquired by the electrodes C3 and C4 respectively have been represented. However, you can get the power spectra just by changing the electrodes according to the 10/20 standard. Although this choice may cause a variation in the pre-processing stage, since it all depends on the nature of the signal.

Motivation:

How to plot an EEG signal?

First, you have load the file called "Subject02.mat" in your version of matlab, then in the Command Window type the following line of code:

>> x_k = Registro(1).Izquierda(:,17,1);
>> N = length(x_k);
>> Fm = 250;
>> T = 1./Fm;
>> n = 0:N-1;
>> k = n*T;
>> plot(k,x_k)

This line graphs the signal captured at electrode C3, first trial and first run. However, it should be noted that this signal presents a trend and random noise, so the signal behaves as a non-stationary time series. To correct this problem, implement the following lines of code:

>> x_k = Registro(1).Izquierda(:,17,1);
>> N = length(x_k);
>> Fm = 250;
>> T = 1./Fm;
>> n = 0:N-1;
>> k = n*T;
>> xd_k = detrend(x_k);
>> FN = Fm/2;
>> wsh = 0.1/FN;
>> wph = 1/FN;
>> Rs = 30;
>> Rp = 5;
>> [Order, frec] = buttord(wph,wsh,Rp,Rs);
>> [b,a] = butter(Order,wph,'high');
>> hph_k = filter(b,a,xd_k);
>> wpl = 55/FN;
>> wsl = 60/FN;
>> [Orderl, frecl] = buttord(wpl,wsl,Rp,Rs);
>> [bl,al] = butter(Orderl,wpl,'low');
>> hpl_k = filter(bl,al,hph_k);
>> figure(1)
>> hold on
>> plot(k,xd_k)
>> plot(k,hpl_k)
>> legend('detrend signal','Hp+Lp filter')
>> hold off

You must note that when applying the signal filtering, a phase shift is obtained which is caused by the frequency response of the filters, that is, if we take an analog circuit as a reference, the capacitive elements of which are make up the firtos cause the signal to have a frequency response, with defined magnitude and phase parameters, that is to say that said signal is transformed into a phasor.

In some applications, the signal delay can become extremely important, which is why this must be corrected. One way to do this is to implement an all-pass filter in the output circuit. In this algorithm the all-pass filter section is not implemented, so we are depressing the delay time of the h [n] signals from the outputs of the alpha and beta bands.

In summary, in the algorithm you will find the pre-processing stage (elimination of trend), conformation of the alpha and beta oscillations in the sensorimotor bands or SensoriMotor Rhytms (SMRs) through the use of elliptical filters and the calculation of the PSD by Welch's method using the relation P (dB) = log (P), all thanks to the recordings made in the electrodes C3 and C4 respectively.

To run the algorithm you must open your version of Matlab, load the file "Sujeto02.mat", load the file called "EEG_PSD.m" and press F5 on your computer keyboard or click on the "Run" button (color green).

Enjoy, consult, learn and advance :)