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.