//// Power spectrum of Izhekevich data
 for i=101:1000; y(i)=sum(firings(:,1)==i); end
 Y = fft(y)/1000;
 N = length(Y); Y(1)=[];
 power = abs(Y(1:N/2)).^2;
 freq = (1:N/2)/N*1000;
 plot(freq,power)
 xlabel('Frequency (Hz)'); ylabel('Power')