%This function cleanses a single cell's waveform, extracts spikes, and
%displays the results.  Suggested variable parameters are given inside
%[ ]'s.
function [wave_cleansed,spks]=DataCleansing_Mult(dataStruct,trial,cell,thresh,BW,LPF_cut,HPF_cut,spk_window,buff_window,end_window)

trial;              %Trial # we wish to look at
cell;               %Cell # we are looking at
thresh;             %Threshold in spectrum to zero out [10]
BW;                 %Bandwidth around peaks to zero out (Hz) [4]
LPF_cut;            %Cutoff of LPF (Hz)...has theoretical justification [5e3]
HPF_cut;            %Cutoff of HPF (Hz)...also has justification [400]
spk_window;         %0's window size for spike done (seconds) [0.4e-3]
buff_window;        %Buffer time before looking for end of spike (seconds) [0.5e-3]
end_window;         %End time to stop looking for end of spike (seconds) [3e-3]
samp_rate = 20000;  %Sampling rate (Hz)
num_dev = 3.2;        %Standard deviation at which spike is detected
chop = 250;         %Index of a time sufficiently past stimulus (used to detect std dev properly)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prep original data & show spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Grab cell of appropriate trial #
wave = dataStruct.data{trial};

%Detrend data -- only care about AC b/c capactively coupled
orig_wav = detrend(wave(:,cell));

%Plot original waveform
figure(1);
subplot(4,1,1);
plot(dataStruct.timeVec{trial},orig_wav,'b-');
ylabel('Amplitude (mV)');
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{trial});
freq_step = 1000/dataStruct.timeVec{trial}(tot_samps); % time vec is in ms!

%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;
xlabel('Freq (Hz)');
title('Spectral Cleansing');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 only AC changes (spikes)
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{trial},wave_cleansed,'r-');
ylabel('Amplitude (mV)');
title('Cleansed Signal');


%%%%%%%%%%%%%%%%%%%
% 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(chop:length(wave_cleansed))));

%Detect large signal deviations (+/- 'num_dev' std_devs)
for i = 1:length(wave_cleansed),
     if(abs(wave_cleansed(i)) >= num_dev*std_dev)
         devs(i) = 1;
     else
         devs(i) = 0;
     end;
end;

%Plot deviations waveform
figure(1);
subplot(4,1,3);
plot(dataStruct.timeVec{trial},devs,'b-');
ylabel('T/F');
msg = sprintf('Deviation > +/- %f sigma',num_dev);
title(msg);

%Allocate memory for spike train
spks = zeros(1,length(wave_cleansed));

% %Determine spike center of mass
% i = 1;
% while(i <= length(devs)),
%      if(devs(i) == 1),
%          mass = 0;
%          sum = 0;
%          for j = i:(i+samp_rate*spk_window),
%              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;

figure(1);
subplot(4,1,4);
plot(dataStruct.timeVec{trial},spks,'g-');
xlabel('Time (ms)');
ylabel('Spike');
title('Extracted Spike Train');

return;