| 
  • If you are citizen of an European Union member nation, you may not use this service unless you are at least 16 years old.

  • You already know Dokkio is an AI-powered assistant to organize & manage your digital files & messages. Very soon, Dokkio will support Outlook as well as One Drive. Check it out today!

View
 

FFT_analysis

Page history last edited by Ki-Young Jung 3 years, 2 months ago

대한뇌파연구회 집담회: 뇌파분석워크샵

 
This tutorial is published online by 
Professor Ki-Young Jung, MD, PhD
ENSLab. Neuroscience Research Institute
Seoul National University College of Medicine
last updated Oct. 18 2019
 
This tutorial is made using a ERP data published in following paper.
Sunwoo JS et al. Abnormal activation of motor cortical network during
phasic REM sleep in idiopathic REM sleep behavior disorder.
Sleep 2019; 42. DOI 10.1093/sleep/zsy227

 

 

Contents

 

make sine wave in 3, 5, 7 Hz

freq = [3 5 7];
srate = 200;
time = -1:1/srate:1;

% 3 Hz sine wave
sine_3Hz = exp(1i*2*pi*freq(1).*time);
figure; plot(time, real(sine_3Hz));

sine_5Hz = exp(1i*2*pi*freq(2).*time);
figure; plot(time, real(sine_5Hz));

sine_7Hz = exp(1i*2*pi*freq(3).*time);
figure; plot(time, real(sine_7Hz));

sum_sine = sine_3Hz + sine_5Hz + sine_7Hz;
figure; plot(time, real(sum_sine));

% add noise
noise = rand(1, length(sum_sine));
sum_sine_noise = sum_sine + noise;
figure; plot(time, real(sum_sine_noise))

%  fourier transform
data = sum_sine;
N = length(data);
nyquist = srate/2;
frequencies = linspace(0,nyquist,N/2+1);
time = ((1:N)-1)/N;

% initialize Fourier output matrix
fourier = zeros(size(data));

% Fourier transform is dot-product between sine wave and data at each frequency
for fi=1:N
    sine_wave   = exp(-1i*2*pi*(fi-1).*time);
    fourier(fi) = sum(sine_wave.*data);
end

fourier=fourier/N;

fft_power = abs(fourier(1:N/2+1)).^2;

% plot

figure;
% plot(frequencies,fft_power, 'r', 'linewidth', 3)
bar(frequencies,fft_power)
set(gca,'xlim',[0 20])
xlabel('Frequency (Hz)')
ylabel('power(uV^2)')

figure;
plot(frequencies,fft_power, 'r', 'linewidth', 3)
set(gca,'xlim',[0 20])
xlabel('Frequency (Hz)')
ylabel('power(dB)')

FFT using REM sleep EEG data

 

Download data

RBD_C12_to.mat

RBD_C12_ph.mat

 

 

1. single epoch 

method 1. using dot product with sine wave

load RBD_C12_to.mat
eegplot(Data); % require EEGLAB to plot raw EEG 

take a single epoch for set parameter

data = Data(18, :, 1);  % ch * data * epoch

calculate FFT

srate = 200;
N = length(data)
nyquist = srate/2;
frequencies = linspace(0,nyquist,N/2+1);
time = ((1:N)-1)/N;

% initialize Fourier output matrix
fourier = zeros(size(data));

% Fourier transform is dot-product between sine wave and data at each frequency
for fi=1:N
    sine_wave   = exp(-1i*2*pi*(fi-1).*time);
    fourier(fi) = sum(sine_wave.*data);
end

fourier=fourier/N;
fft_power = abs(fourier(1:N/2+1)).^2;

% plot PSD
figure;
plot(frequencies,fft_power, 'r', 'linewidth', 2)
set(gca,'xlim',[0 30])
xlabel('Frequency (Hz)')
ylabel('power(uv^2)')

% band power
alpha_idx = dsearchn(frequencies', [8, 13]')
alpha_power = sum(fft_power(alpha_idx(1):alpha_idx(2)))

method 2. using fft function

fft_data = fft(data);
fft_data = fft_data / N;

fft_power = abs(fft_data(1:N/2+1)).^2;

figure;
plot(frequencies, fft_power, 'r', 'linewidth', 2)
set(gca,'xlim',[0 30])
xlabel('Frequency (Hz)')
ylabel('power(uV^2)')

2. whole epoch: phasic REM sleep

load data and set parameter

load RBD_C12_ph.mat
epoch = Data(18, :, 1);  % ch * data * epoch
srate = 200;
N = length(epoch);
nyquist = srate/2;
frequencies = linspace(0,nyquist,N/2+1);

% calculate fft
for epoch =1:size(Data, 3)

data = Data(5, :, epoch);  % ch * data * epoch
fft_data = fft(data);
fft_data = fft_data / N;
fft_power(epoch, :) = abs(fft_data(1:N/2+1)).^2;

end

mean_fft_power_ph = mean(fft_power);

figure;
plot(frequencies, mean_fft_power_ph, 'r', 'linewidth', 2)
set(gca,'xlim',[0 30])
xlabel('Frequency (Hz)')
ylabel('power(uV^2)')

whole epoch: tonic REM sleep %%

load RBD_C12_to.mat

epoch = Data(18, :, 1);  % ch * data * epoch
srate = 200;
N = length(epoch);
nyquist = srate/2;
frequencies = linspace(0,nyquist,N/2+1);

for epoch =1:size(Data, 3)

data = Data(5, :, epoch);  % ch * data * epoch
fft_data = fft(data);
fft_data = fft_data / N;
fft_power(epoch, :) = abs(fft_data(1:N/2+1)).^2;

end

mean_fft_power_to = mean(fft_power);

figure;
plot(frequencies, mean_fft_power_to, 'b', 'linewidth', 2)
set(gca,'xlim',[0 30])
xlabel('Frequency (Hz)')
ylabel('power(uV^2)')

compare phasic with tonic REM sleep

figure;  hold on
plot(frequencies, mean_fft_power_ph, 'r', 'linewidth', 2)
set(gca,'xlim',[0 30])
xlabel('Frequency (Hz)')
ylabel('power(uV^2)')

plot(frequencies, mean_fft_power_to, 'b', 'linewidth', 2)
set(gca,'xlim',[0 30])
xlabel('Frequency (Hz)')
ylabel('power(uV^2)')

legend('phasic', 'tonic')

psd for averaged eeg data

avg_data = mean(data, 3);
fft_data = fft(avg_data);
fft_data = fft_data / N;
fft_power = abs(fft_data(1:N/2+1)).^2;

figure;
plot(frequencies, fft_power, 'k', 'linewidth', 2)
set(gca,'xlim',[0 30])
xlabel('Frequency (Hz)')
ylabel('power(uV^2)')

Comments (0)

You don't have permission to comment on this page.