% This code assumes openEEG.m has already been run. Try it with N set to larger
% values to get a good set of data (e.g. N = 256*30 gives 30 seconds)
% Tim Collins, 2010

% First perform time-frequency analysis using STFT
% This example uses data from channel 1 and a window size of 256 samples (1
% second). Longer windows will give better frequency resolution but worst
% time resolution.
channel = 1;
windowSize = 256;

%[y,f,t,p] = spectrogram(eegData(:,channel),windowSize,[],[],256);
%Work out spectrogram without using signal processing toolbox...
M = 7680 * 2 / windowSize - 1;
x = zeros(windowSize, M);
win = 0.54 - 0.46 * cos(2*pi*(0:(windowSize-1))'/(windowSize-1));  % Hamming window
for n = 1:M
    x(:,n) = win .* eegData((1:windowSize) + (n-1)*windowSize/2, channel);
end
y = fft(x);
p = abs(y(1:(windowSize/2 + 1),:) / sum(win)).^2 *sqrt(2);
p(1,:) = p(1,:) * .5;
p((windowSize/2 + 1),:) = p((windowSize/2 + 1),:) * .5;
f = (0:(windowSize/2))';
t = (0:(M-1)) * windowSize / 512;

% Now plot the results. I've restricted the frequency range from 0-40 Hz
% here.
surf(t,f(1:41),10*log10(abs(p(1:41,:))),'EdgeColor','none');
axis xy; axis tight; colormap(jet); view(0,90);
xlabel('Time');
ylabel('Frequency (Hz)');
title('Power spectral density [dB]');
colorbar vert

Published with MATLAB® 7.4