%This function cleanses all cells from single trial data and combines the
%results into a single file.
function SpikeExtract(file_name)

thresh = 10;            %Threshold in spectrum to zero out
BW = 4;                 %Bandwidth around spectral peaks to zero out (Hz)
LPF_cut = 5e3;          %Cutoff of LPF (Hz)...has theoretical justification
HPF_cut = 400;          %Cutoff of HPF (Hz)...also has justification
spk_window = 0.4e-3;    %0's window size for spike done detection (seconds)
buff_window = 0.5e-3;   %Buffer time before looking for end of spike (seconds)
end_window = 3e-3;      %End time to stop looking for end of spike (seconds)
samp_rate = 20000;      %Sampling rate (Hz)
num_devs = 3.2;           %Look for time-domain excursions exceeding this # std_devs

pts_per_file = zeros(6,1);      %Record # samples per file, for later

%For the six broken-up files
for file_num = 1:6,
        
	%Load file
	file = sprintf('%s%d.mat',file_name,file_num);
	load(file);
    
    pts_per_file(file_num) = length(dataStruct.timeVec);
    
    %Pre-allocate memory
	spk_trains = zeros(size(dataStruct.data));
    
    %Notify user
    msg = sprintf('Extracting data set %d...',file_num);
    disp(msg);
	
        %For each cell
		for cell = 1:64,
            
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		% Prep original data
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
		%Detrend data -- only care about AC b/c capactively coupled
		orig_wav = detrend(dataStruct.data(:,cell));
		
		% %Plot original waveform
		% figure(1);
		% subplot(4,1,1);
		% plot(dataStruct.timeVec,orig_wav,'b-');
		% xlabel('Time (ms)');
		% ylabel('Amplitude');
		% title('Original signal');
		
		%Take FFT
		spec = fft(orig_wav);
		
		%Calculate frequency step in FFT (df = 1/total time length of data)
		%One step in FFT array = df in frequency
		tot_samps = length(dataStruct.timeVec);
		freq_step = 1000/dataStruct.timeVec(tot_samps); % time vec is in ms!
		time_vec = dataStruct.timeVec;
        
 		% %Construct a freq vector corresponding to 1/2 of FFT data (other 1/2 mirrors)
 		% freq = [0:freq_step:tot_samps*freq_step/2];
 		
		% %Plot one-sided spectrum of original data (abs b/c we only want magnitude)
		% figure(2);
		% plot(freq,abs(spec(1:(tot_samps/2)+1)),'b-');
		% hold on;
		
		
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		% Adaptive spectral zeroing of deterministic signals
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
		%Local variables
		f_zero = [];        %Spectral loctaions to zero out
		
		
		%Find spectral peaks above thresh and add to f_zero
		%Note: takes care of two-sided also, b/c both
		%mirrored peaks are flagged
		for i = 1:length(spec),
            if(abs(spec(i)) > thresh)
                f_zero = [f_zero,i];
            end;
		end;
		
		%Zero out appropriate locations in the spectrum
		%Convert BW into approprite number of points in FFT!!
		for i = 1:length(f_zero),
            
            %Calculate window to be zeroed
            fl = round(f_zero(i)-BW/2/freq_step);
            fh = round(f_zero(i)+BW/2/freq_step);
            
            %If window exceeds available freqs, correct it
            if (fl < 1),
                fl = 1;
            elseif(fh > length(spec)),
                fh = length(spec);
            end;        
            
            %Null out frequency components in FFT
            spec(fl:fh)=0;
		
		end;
		
		%LPF signal to remove extraneous spectral data
		i_low = round((LPF_cut/freq_step)+1);
		i_high = round(length(spec)-LPF_cut/freq_step+1);
		spec(i_low:i_high) = 0;
		
        
		%HPF signal to increase response to AC changes (spikes)
        %and decrease to low-frequency fluctuations
		i_low = round((HPF_cut/freq_step)+1);
		i_high = round(length(spec)-HPF_cut/freq_step+1);
		spec(1:i_low)=0;
		spec(i_high:length(spec))=0;

		% %Superimpose new spectrum over old...peaks are gone & LPF
		% plot(freq,abs(spec(1:(tot_samps/2)+1)),'r-');
		
		%Reconstruct filtered waveform
		%Need only real b/c imag is zero...as expected
		wave_cleansed = real(ifft(spec));
		wave_cleansed = wave_cleansed - mean(wave_cleansed);
		
		% %Plot cleaned waveform
		% figure(1);
		% subplot(4,1,2);
		% plot(dataStruct.timeVec,wave_cleansed,'r-');
		
		
		%%%%%%%%%%%%%%%%%%%
		% Extract spikes
		%%%%%%%%%%%%%%%%%%%
		
		%Must preallocate to avoid memory issues!
		devs = logical(zeros(1,length(wave_cleansed)));
		
		%Find std. dev. of waveform
		std_dev = sqrt(var(wave_cleansed));
		
		%Detect large signal deviations (+/- 'num_devs' std_devs)
		for i = 1:length(wave_cleansed),
             if(abs(wave_cleansed(i)) >= num_devs*std_dev)
                 devs(i) = 1;
             else
                 devs(i) = 0;
             end;
		end;
		
		%Allocate memory for spike train
		spks = zeros(length(wave_cleansed),1);
		
% 		%Determine spike center of mass
% 		i = 1;
% 		while(i <= length(devs)),
%              if(devs(i) == 1),
%                  mass = 0;
%                  sum = 0;
%                  for j = i:min((i+samp_rate*spk_window),length(devs)),
%                      if(devs(j) == 1),
%                         mass = mass + j;
%                         sum = sum + 1;
%                      end;
%                  end;
%                  spks(floor(mass/sum)) = 1;
%                  i = i+samp_rate*spk_window+1;
%              else
%                  i = i + 1;
%              end;
% 		end;

        
        %Determine spike locations, spike over if 0 for > spk_window
		i = 1;
		while(i <= length(devs)),
             if(devs(i) == 1),
                %Found starting edge of spike
                spks(i) = 1;
                 
                %Search for end of spike
                %Give buff_window starting buffer, and search until end_window later
                j = i + samp_rate*buff_window;
                consec_low = 0;
                while((j <= (i+samp_rate*end_window)) && (j <= length(devs)) && (consec_low <= spk_window*samp_rate)),
                    if(devs(j) == 0),
                        consec_low = consec_low + 1;
                        j = j + 1;
                    else
                        consec_low = 0;
                        j = j + 1;
                    end;
                end;
                
                %Set i for continued spike detection
                i = j;
            else
                %Else look at next location
                i = i + 1;
            end;
		end;
        
        %Store spks in primary matrix	
		spk_trains(:,cell) = spks;
		
end;
	
    %Convert spike train matrix to sparse matrix
    [i,j,s] = find(spk_trains);
    [m,n] = size(spk_trains);
    spk_trains = sparse(i,j,s,m,n);
    
    %Save extracted spike trains to file
	file = sprintf('%s%d_spks.mat',file_name,file_num);
	save(file,'spk_trains','time_vec');
    
    %Temporarily deallocate memory
    clear spk_trains dataStruct;
    
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Combine data into one file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Compressing into 1 file...');

%Allocate memory
tot_pts = 0;
for i = 1:6,
    tot_pts = tot_pts + pts_per_file(i);
end;
extracted_spks = logical(zeros(tot_pts,64));
time = zeros(tot_pts,1);

%Determine splice points
splices = zeros(7,1);
tmp = 0;
for i = 2:7,
    tmp = tmp + pts_per_file(i-1);
    splices(i) = tmp; 
end;

%Pull in spike data from disk and add to matrix
for file_num = 1:6,
    
    %Load data
    file = sprintf('%s%d_spks.mat',file_name,file_num);
	load(file);
    %Update matrices
    extracted_spks((splices(file_num)+1):splices(file_num+1),:) = spk_trains(:,:);
    time((splices(file_num)+1):splices(file_num+1),1) = time_vec(:,:);
    
end;

%Convert spike train into sparse matrix
[i,j,s] = find(extracted_spks);
[m,n] = size(extracted_spks);
extracted_spks = sparse(i,j,s,m,n);

%Save single spike train file
file = sprintf('%s_extracted.mat',file_name);
save(file,'extracted_spks','time');
    
return;