Matlab - Espectro de señal

   
Vista:

Espectro de señal

Publicado por ivan (5 intervenciones) el 27/08/2014 00:35:58
Hola que tal, soy nuevo por esta comunidad...
Estoy trabajando con señales EEG (electroencefalograma), a dichas señales les aplico el filtro FIR pasa altas y posteriormente la FFT pero creo que me esta fallando la FFT ya que la gráfica para todos los canales que tengo (14) siempre es la misa.

Alguna sugerencia de como hacerlo correctamente?

Como puedo obtener el espectro para ondas alpha (8 - 13 Hz) con la FFT?


Este es el código que estoy utilizando
muestras por segundo: 128
cada 0.0078125 s se toma una muestra

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
%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


Saludos y muchas gracias por la ayuda.

untitled
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder

Espectro de señal

Publicado por khin Rosales (1 intervención) el 15/11/2016 02:53:49
Hola buen dia, quisa yo te pueda ayudar solo dame un poco de tiempo y si me permites usar tu codigo con una funcion para una clase de procesamiento digital de senales y si lo corrigo lo subire aqui mismo.

Saludos.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por Joshua (4 intervenciones) el 01/12/2016 07:06:37
Buenas, estoy trabajando en un caso muy parecio al tuyo, pero tengo problemas al momento de hacer el espectrograma de mi señal, quisiera saber si podiste solucionarlo, te lo agradeceria mucho.
Atte: Joshua
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por Joshua (4 intervenciones) el 01/12/2016 07:07:08
Buenas, estoy trabajando en un caso muy parecio al tuyo, pero tengo problemas al momento de hacer el espectrograma de mi señal, quisiera saber si podiste solucionarlo, te lo agradeceria mucho.
Atte: Joshua
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por Ivan (5 intervenciones) el 01/12/2016 07:19:02
Hola yoshua, claro, no hay problema
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por Ivan (5 intervenciones) el 01/12/2016 07:24:29
@yoshua, si pude solucionar el problema que tenía para que me saliera de manera adecuada, de echo el código es diferente al que indique aquí, en una chanza que tenga publico la manera en lo que solucione
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por Joshua (4 intervenciones) el 01/12/2016 07:34:01
Gracias Iván , te lo agradecería mucho es para un proyecto el cual estoy trabajando o si pudieras mandarmelo a mi correo te lo agradecería muchísimo joshua.aragon@ucsp.edu.pe
Saludos
Atte: Joshua
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por ivan (5 intervenciones) el 05/12/2016 05:47:12
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
@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');
%}
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por Joshua (4 intervenciones) el 05/12/2016 07:48:51
Hola Ivan, muchas gracias por tu ayuda!!! pero disculpa que tengo unas consultas, que no me quedaron claro y te agradeceria si me las pudieras aclarar.
Pd: disculpa si las preguntas son algo bobas, pero gracias por tu tiempo.

No me queda muy claro esta parte:
1
2
3
4
5
6
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;

Un poco mas arriba al sacar la fft, porque amplias su tamaño de la señal a N??
Mas abajo, con que fin sacas la magnitud y la potencia

Muchas gracias por tu tiempo, saludos.
Atte: Joshua
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Espectro de señal

Publicado por Ivan (5 intervenciones) el 05/12/2016 19:46:59
@Joshua, hola que tal, en las lineas
1
2
3
4
5
6
7
8
9
10
11
12
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;

lo que hago es aplicar un ventaneo a la señal, checa cuales hay (no recuerdo el nombre ahorita) pero hay como 5 muy utilizados, cada uno tiene una formula.
La magnitud y la potencia para hacer comparaciones entre diferentes segmentos para ver si hay señales parecidas (en magnitud y potencia) y para graficar cada tramo utilizado.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar