EEG DWT : How to get desired frequency spectrum after DWT

5 views (last 30 days)
Hi, I am working on EEG signal. At first i applied the Butterworth Low Pass Filter to extract 0-64 Hz frequency. Then i applied DWT to extract BETA (16-32Hz) and ALPHA(8-16Hz) wave . So, according to theory , 2nd and 3rd level coefficient of DWT should provide the beta and alpha wave. But when i performed the fft of D3 wave i did not get the desired spectrum. I cant understand the problem. Please help.......
Following is my code. I used the EEG DATASET-1 of the following link....... https://sites.google.com/site/projectbci/
MATLAB CODE:
clc
clear all
load('C:\Users\COOL\Documents\MATLAB\Subject1_1D.mat');
fs=500;
T=1/fs; %T=.002
tmax=3;
Nsamps = tmax*fs;
t = 0:1/fs:tmax-1/fs;
s = right(6,1:1500);
l_s = length(s);
%Plot in Time Domain Original EEG
figure(1)
plot(t,s)
xlabel('Time (s)')
ylabel('Amplitude (V)')
title('Original Signal')
ylim([-100 100])
grid
%Spectrum of origonal EEG
N = fs;
S=fft(s,N);
CS=[S(1:N/2+1)];
freq=[0:N/2]/(N*T);
mag=abs(CS);
figure(2);
plot(freq,mag);
title('Spectrum of Original Signal')
xlabel('Frequency (Hz)');
ylabel('Amplitude');
ylim([0 1000])
grid
%Apply butterworth Low Pass Filter
N=6;
Wn_t = 64/fs;
[c,d] = butter(N,Wn_t);
s1 = filtfilt(c, d, s);
%Spectrum of filtered EEG N = fs;
S=fft(s1,N);
CS=[S(1:N/2+1)];
freq=[0:N/2]/(N*T);
mag=abs(CS);
figure(3);
plot(freq,mag);
title('Spectrum of filtered Signal')
xlabel('Frequency (Hz)');
ylabel('Amplitude');
ylim([0 1000])
grid %Perform DWT
[C,L] = wavedec(s1,3,'db4');
cA5 = appcoef(C,L,'db4',3);
cD3 = detcoef(C,L,3);
cD2 = detcoef(C,L,2);
cD1 = detcoef(C,L,1);
A3 = wrcoef('a', C, L,'db4',3);
D3 = wrcoef('d',C,L,'db4',3);
D2 = wrcoef('d',C,L,'db4',2);
D1 = wrcoef('d',C,L,'db4',1);
%Spectrum of filtered D3(Level 3 coefficient reconstructed)
N = fs;
S=fft(D3,N);
CS=[S(1:N/2+1)];
freq=[0:N/2]/(N*T);
mag=abs(CS);
figure(4);
plot(freq,mag);
title('Spectrum of D3')
xlabel('Frequency (Hz)');
ylabel('Amplitude');
ylim([0 500])
xlim([0 100])
grid

Answers (0)

Categories

Find more on EEG/MEG/ECoG in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!