FindVPAamp(vpa)

function Processed = FindVPAamp(vpa)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function identifies local max and min for amplitude calculation and
% averages in times bins. Requested time bins are calculated by 
% user-supplied sampling rate (Hz). Uses zero-crossing assumption 
% following detrend. Artifacts should be minimized PRIOR to using this algorithm
%
% Input:    Matrix with subjects (or condition/subject) on rows
%           Sampling rate (Hz) manually
%           Desired bins of time (s) manually
% Output:   .dat file with average values per bin per row
%           .dat file with every peak amplitude (allpeaks.dat)
%
% Developement goals:
% - Add check that file is as long as necessary (or cut) for bin length
% - Add figures for allpeaks to visually identify processing problems
% Author: Nicole Prause 12.15.08
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%For wavelet project looping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% delrows=2:1:43; %As matrix collapses will not actually need to delete "every other"
% for gottago = 1:length(delrows)
%     signal(delrows(gottago),:)=[];
% end

format short G

output=[];
Bin = input('Enter how many seconds you want in each bin:'); %Time length of bins requested
SRate = input('Enter the sampling rate (in Hz) of your signal:'); %Sampling rate
SRate=80;

[r,c] = size(vpa); %Get number of rows to set up how many to run

for row = 1:r
   signal=vpa(row,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Record time of zero-crossings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signal = detrend(signal); %Removes linear drift; zero-crossing dependent
start=signal(1:end-1);
stop=signal(2:end);
time=start.*stop; %Compares if sign of closest available (restricted by sampling rate) data points change
cross=find(time<0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate amplitude per every two (peak/trough) zero-crossings
% Uses non-overlapping windows
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:2:(length(cross)-2)
    StartPeak=cross(i);
    Middle=cross(i+1);
    EndPeak=cross(i+2);
    PeakBin=signal(1,StartPeak:Middle);
    TroughBin=signal(1,Middle:EndPeak);
    AmplitudePeak=max(abs(PeakBin));
    AmplitudeTrough=max(abs(TroughBin));
    Amplitude=AmplitudePeak + AmplitudeTrough;
    AllPeakAmplitudes=[StartPeak Amplitude];

    if i==1
        dlmwrite('allpeaks.dat', AllPeakAmplitudes)
    else
        dlmwrite('allpeaks.dat', AllPeakAmplitudes, '-append')
    end
end
    load allpeaks.dat
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Average peaks in bins
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CalcBin = Bin*SRate; %Length (in data points) per bin
BinMark = 1:CalcBin:length(signal); %Data point for start of each bin

for j=2:length(BinMark)
        l=find(allpeaks(:,1)>BinMark(j-1) & allpeaks(:,1)<BinMark(j));

AvgVector = [];
    for index=1:length(l)
        AvgVector(index,:)= allpeaks(l(index),2);
    end

    Avg= mean(AvgVector);
    
    if j==2
        dlmwrite('averages.dat', Avg)
    else
        dlmwrite('averages.dat', Avg, '-append')
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load the processed data for exporting, processing as-needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load averages.dat
averages = averages';

if j==2
        dlmwrite('output.dat', averages)
    else
        dlmwrite('output.dat', averages, '-append')
end

load output.dat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Option to inspect averages
% - Requires manual intervention to advance
% - May be removed without consequence to processing routine
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(output(row,:));
title(['Averages for subject (or condition) in row ',num2str(row)])
xlabel ('Time bin')
ylabel ('VPA average')
display('Press any key to continue')
pause

end