Matlab - Espectro de señal

%se carga el archivo csv y se asignan los datos a data

data = csvread('csv.csv', renglon, 0);%csvread('csv.csv', row, col)

%se asignana los datos a cada canal

tiempo = data(:,1);

AF3 = data(:,3);

F7 = data(:,4);

F3 = data(:,5);

FC5 = data(:,6);

T7 = data(:,7);

P7 = data(:,8);

O1 = data(:,9);

O2 = data(:,10);

P8 = data(:,11);

T8 = data(:,12);

FC6 = data(:,13);

F4 = data(:,14);

F8 = data(:,15);

AF4 = data(:,16);

cont = 130;

%cambair a secuencia los numeros de muestra

for t = (cont-1):length(tiempo)-1

    tiempo(cont) = t;

    cont = cont+1;

end

%multiplicar cada muestra por 0.0078125 (tiempo en que se toma cada muestra)

for t = 1:length(tiempo)

    tiempo(t) = (tiempo(t)*0.0078125);

end

%APLICANDO FILTRO FIR PASA AILTAS

sample_rate = 128; %frecuencia de muestreo

nsamples = 512; %muestras

t1 =([tiempo(1:512)]); %t = ([F4(:512)]-1);

% Choose filter cutoff frequency ()

cutoff_hz = 0.2;

% Normalize cutoff frequency (wrt Nyquist frequency)

nyq_freq = sample_rate / 2; %calcular el numero de Nyquist

cutoff_norm = cutoff_hz / nyq_freq;

% FIR filter order (i.e. number of coefficients - 1)

order = 2;

fir_coeff = fir1(order, cutoff_norm,'high');

% Filter the signal with the FIR filter

filtered_signal = filter(fir_coeff, 1, F3);

plot(tiempo(1:128),filtered_signal(1:128),'b',tiempo(1:128),O1(1:128),'r')

%--------------

% Group delay as a scalar, in number of samples

group_delay = median(grpdelay(fir_coeff));

% Group delay in seconds

group_delay_s = group_delay / sample_rate;

% Plot the original signal...

figure(1)

plot(t1(1:512), O1(1:512))

% ...and allow adding more plots

hold on

% Align and plot the filtered signal

% (On top of existing one)

plot(t1(order:end)-group_delay_s, filtered_signal(order:512), 'g')

% Align and plot only the d6esired part of the filtered signal (discarding

% the transient)

plot(t1(order:end)-group_delay_s, filtered_signal(order:512), ...  %plot(t1(order:end)-group_delay_s, filtered_signal(order:end), ...

    'r');

grid on

hold on

%ESTA PARTE ES DE LA FFT

%serie temporal

n=2^14; % 2^12 = 4096 datos

dt=0.0078125; %intervalo de tienpo (muestreo)

t=(0:n-1)*dt; %vector de tiempos (0:n-1)*dt

%amplitud-fase vs. frecuencias

g=fft(filtered_signal); %g=fft(x);

power=abs(g).^2;

dw=2*pi/(n*dt);

w=(0:n-1)*dw; %vector de frecuencias angulares

plot(w,power(1:length(w)),'m')

xlabel('\omega')

ylabel('P(\omega)')

title('Espectro de potencia')

grid on

hold off

@Joshua, espero te sirba de algo, saludos

%data = eegdata;                %señal EEG

%se separa la señal por electrodos

AF3 = data(:,1);

F7 = data(:,2);

F3 = data(:,3);

FC5 = data(:,4);

O1 = data(:,5);

O2 = data(:,6);

FC6 = data(:,7);

F4 = data(:,8);

F8 = data(:,9);

AF4 = data(:,10);

fs = 128;                      %frecuencia de muestreo

N = 2^nextpow2(length(data));  %Next power of 2 from length of data

n = length(data);              %tamaño de la señal (datos)

%se remueve la media de la señal por cada electrodo

AF3 = AF3 - mean(AF3);

F7  = F7 - mean(F7);

F3  = F3 - mean(F3);

FC5 = FC5 - mean(FC5);

O1  = O1 - mean(O1);

O2  = O2 - mean(O2);

FC6 = FC6 - mean(FC6);

F4  = F4 - mean(F4);

F8  = F8 - mean(F8);

AF4 = AF4 - mean(AF4);

data = [AF3, F7, F3, FC5, O1, O2, FC6, F4, F8, AF4];  %se agrupan los electrodos a la matriz

eeg = detrend(data);           %remueve los valores medios / tendencia  lineal

%Se crea el filtro Butterworth

[a, b] = butter(6, 29.86/fs);        %corte ssuperior del filtro 35Hz

[c, d] = butter(6, 0.85/fs, 'high'); %corte inferior del filtro 1Hz

%se aplica el filtrado a la señal

eeg = filter(a, b, eeg);      %se aplica el corte superior

eeg = filter(c, d, eeg);      %se aplica el corte inferior

pi = 3.1415;

npts = length(N);

a = npts/128;

x = linspace(0,a,npts);

A = 0.54-0.46*(cos(2*pi(x/a)));

eeg = A'*eeg;

%se aplica la FFT a la señal filtrada

eegFFT = fft(eeg, N);

%guardar la señal procesada

csvwrite('eeg.csv',eegFFT);

%se obtiene el poder espectral (parte real de los datos FFT)

%separo los valores por electrodo

AF3 = eegFFT(:,1);

F7 = eegFFT(:,2);

F3 = eegFFT(:,3);

FC5 = eegFFT(:,4);

O1 = eegFFT(:,5);

O2 = eegFFT(:,6);

FC6 = eegFFT(:,7);

F4 = eegFFT(:,8);

F8 = eegFFT(:,9);

AF4 = eegFFT(:,10);

%obtengo la magnitud de cada electrodo

AF3 = abs(AF3);

F7 = abs(F7);

F3 = abs(F3);

FC5 = abs(FC5);

O1 = abs(O1);

O2 = abs(O2);

FC6 = abs(FC6);

F4 = abs(F4);

F8 = abs(F8);

AF4 = abs(AF4);

eegH = [AF3, F7, F3, FC5, O1, O2, FC6, F4, F8, AF4];

%guardar la magnitud señal procesada

csvwrite('magnitud.csv',eegH);

%obtengo la potencia de cada electrodo

AF3 = (AF3).^2;

F7 = (F7).^2;

F3 = (F3).^2;

FC5 = (FC5).^2;

O1 = (O1).^2;

O2 = (O2).^2;

FC6 = (FC6).^2;

F4 = (F4).^2;

F8 = (F8).^2;

AF4 = (AF4).^2;

%eegProces = [AF3 F7 F3 FC5 O1 O2 FC6 F4 F8 AF4];

power = [AF3 F7 F3 FC5 O1 O2 FC6 F4 F8 AF4];

%guardar power la señal procesada

csvwrite('power.csv',power);

%power = abs(eegFFT(1:N/2,:)).^2;

%power = abs(eegFFT(1:N/2)).^2;

frecuencia = 0:fs/N:(N/2-1)*fs/N;

%guardar las frecuecnias

csvwrite('frecuencias.csv',frecuencia);

%{

%vector donde se encuentran todas las frecuencias

vectorFrecuencias = [0, 0, ... 

                    ];

for a=1:1:length(power)-1

    vectorFrecuencias(a) = (a)*128/4096;

end

%}

%{

%graficando el resultado

plot(frecuencia,power);

xlabel('Frecuencia Hz');

ylabel('Power');

title('Señal EEG');

%}

pi = 3.1415;

npts = length(N);

a = npts/128;

x = linspace(0,a,npts);

A = 0.54-0.46*(cos(2*pi(x/a)));

eeg = A'*eeg;

pi = 3.1415;

npts = length(N);

a = npts/128;

x = linspace(0,a,npts);

A = 0.54-0.46*(cos(2*pi(x/a)));

eeg = A'*eeg;