PCA filter on EEG data

brain.m

#/usr/bin/octave
%
% Map environment from brain scans
% 
% http://headit-beta.ucsd.edu/recording_sessions/3a7b9ff6-385a-11e3-b29d-0050563f2612
%
% PCA filter the first 2 seconds.
%
% Eigen-decomposing the covariance matrix is the least efficient but easiest way to do PCA.

load 'brain1.mat'

ndata = EEG.data(:, 1:1024)/max(max(EEG.data));

avg = mean(ndata, 1)';
variance = cov(ndata);

[v, lambda] = eig(variance);

% f = ( (xi-m)'C )^2
%
% df/dC = 2 C (xi-m) * (xi-m)  ( (xi-m)'C )
%
% 

%plot(log(abs(lambda*ones(1024,1))));

%pause(5);

a = (0:1:1023)';
% Guassian window in time domain
w = exp(-0.5*(a-512).^2/(2*(1024/12)^2));

% High pass filter in frequency domain
lpf_f = (1-exp(-a));

% Low pass filter in frequency domain
hpf_f = (ones(1024,1) - ones(1024, 1)./(ones(1024, 1)+exp(0.2*(-a+(30/EEG.srate)*1024))));

adjust = eye(1024);
vnew = eye(1024);

for i=1:1024
  printf('%f\n', lambda(end-i+1, end-i+1));
  vi_fft = fft(v(:, end-i+1));
  vnewcol = abs(ifft(vi_fft.*hpf_f.*lpf_f));
  for j = 1:1024
    vnew(j, i) = vnewcol(j);
  end

  if lambda(i, i) <= 1e-3
%    for j = i:1024
%      adjust(end-j+1, end-j+1) = 0;
%    end
%    break;
    adjust(end-i+1, end-i+1) = 0;
  else
    adjust(end-i+1, end-i+1) = 1;
  end
end

newavg = abs(ifft(fft(avg).*hpf_f.*lpf_f));

printf('Done\n');

xr = vnew*adjust*(ndata'-(avg*ones(1, 32)))
xnew = xr %+ (newavg*ones(1, 32));

size(xnew)

for i = 1:32
  %plot(EEG.times(1:1024)*512/1024, abs(fft(xnew(:, i))));
  plot(xnew(:, i));
  pause(2);
end