Comparative_MEA_dataset/Codes/MATLAB spike detection/BatchAnalyser_1.2/FileInput/Axion MATLAB Files/BlockVectorLegacyLoader.m

classdef BlockVectorLegacyLoader < BlockVectorData
    %BLOCKVECTORLEGACYLOADER Static class with resources to load file data
    %   to legacy style structures.
    
    properties(GetAccess = private, Constant = true)
        % Fixed by file format - there's only 1 byte for each
        MAX_CHANNEL_ARTICHOKE = 256;
        MAX_CHANNEL_INDEX     = 256;
    end
    
    methods (Static = true)
        function aData = Legacy_Load_Raw_v1(aSourceData, aHeader, aTimeRange)
            
            % Check to make sure the file format makes sense for a raw file
            if aHeader.header.numSamplesPerBlock ~= 1
                error('load_AxIS_raw:argNumSamplesPerBlock', ...
                    'Invalid header for RAW file: incorrect samples per block');
            end
            
            if aHeader.header.blockHeaderSize ~= 0
                error('load_AxIS_raw:argBlockHeaderSize', ...
                    'Invalid header for RAW file: incorrect block header size');
            end
            
            fseek(aSourceData.FileID, int64(aSourceData.Start), 'bof');
            
            % Copy input structure to output; then add information to output
            aData = aHeader;
            aData.numChannels = aData.header.numChannelsPerBlock;
            if aData.numChannels ~= aData.channelArray.numChannels
                error('load_AxIS_raw:mismatchNumChannels', ...
                    'Invalid RAW file: mismatched number of channels');
            end
            
            if isempty(aHeader.loadedChannels)
                % No channels requested - don't read any data
                return;
            end
            
            % Read data
            fSampleFreq =  int64(aData.samplingFrequency);
            fChannelCount = int64(aData.numChannels);
            fBytesPerSecond = fSampleFreq * fChannelCount * int64(2);
            fMaxTime = int64(aSourceData.EntryRecord.Length / fBytesPerSecond);
            if ~strcmp(aTimeRange, 'all')
                fStart = aTimeRange(1);
				fEnd =aTimeRange(2);
                
                if(fStart >= fEnd)
                    warning('Invalid timespan argument: end time < start time. No valid waveform can be returned.');
                    return;
                end
                if(fStart > fMaxTime)
                    warning( 'DataSet only contains %d Seconds of data (%d Seconds was the requested start time). No waveform will be returned.', fMaxTime, fStart);
                    return;
                end
                if(fEnd > fMaxTime)
                   fEnd = fMaxTime; 
                   warning('DataSet only contains %d Seconds of data. Returned waveforms will be shorter than requested (%d Seconds).', fMaxTime, fEnd - fStart);
                end
                
                fSkipInitialSamples = fStart * fSampleFreq; 
                fSkipInitialBytes = fSkipInitialSamples * fChannelCount * int64(2);
                
                fNumSamples  = int64((fEnd - fStart) * fSampleFreq);
                
				fseek(aSourceData.FileID, fSkipInitialBytes, 'cof');				
            else
                fNumSamples = int64(fMaxTime * fSampleFreq);
            end
            
            if length(aHeader.loadedChannels) == 1
                % We're only reading one channel. For efficiency, take advantage of fread's
                % 'skip' argument.
                
                % First, read and throw away other channels at the beginning of the data
                fread(aSourceData.FileID, aHeader.loadedChannels - 1, 'int16=>int16');
                
                % Read the data for the given channel, skipping from one to the next
                aData.channelData = fread(aSourceData.FileID, fNumSamples, 'int16=>int16', 2*(aData.numChannels - 1));
            else
                % We're reading multiple channels.  There's no easy way to translate this into a
                % single fread, so we read all of the channels in the sample range in question,
                % then remap, throwing away what we didn't want.
                
                % First, read everything within the time range
                fNumSamples      = fNumSamples * int64(aData.numChannels);
                fTempChannelData = fread(aSourceData.FileID, fNumSamples, 'int16=>int16');
                
                % This test (which can fail only when we read the the end of the file)
                % makes sure that we didn't get a number of samples that's not divisible
                % by the number of channels.
                fRemainderCount = mod(length(fTempChannelData), aData.numChannels);
                if fRemainderCount ~= 0
                    warning('load_AXiS_raw:remainderCheck', ...
                        'File %s has wrong number of samples for %u channels', ...
                        aSourceData.FileID, aData.numChannels);
                end
                
                % Convert the 1D array to a 2D array, with channel as the second dimension
                % (starts as first dimension and then is transposed)
                fTempChannelData = reshape(fTempChannelData(:), aData.numChannels, []);
                
                % Remap
                aData.channelData = fTempChannelData(aData.loadedChannels, :);
                aData.channelData = aData.channelData';
            end
        end
        
        function aData = Legacy_Load_spike_v1(aSourceData, aHeader, aTimeRange)
            % Check to make sure the file format makes sense for a spike file
            if aHeader.header.numChannelsPerBlock ~= 1
                error('Load_spike_v1:argNumChannelsPerBlock', ...
                    'Invalid header for SPIKE file: incorrect channels per block');
            end
            
            if aHeader.header.blockHeaderSize < Spike_v1.LOADED_HEADER_SIZE
                error('Load_spike_v1:argBlockHeaderSize', ...
                    'Invalid header for SPIKE file: block header size too small');
            end
            
            if aHeader.header.numSamplesPerBlock < 1
                error('load_AxIS_spike:argNumSamplesPerBlock', ...
                    'Invalid header for SPIKE file: number of samples per block < 1');
            end
            
            % Copy input structure to output; then add information to output
            aData = aHeader;
            
            if isempty(aHeader.loadedChannels)
                % No channels requested - don't read any data
                return;
            end
            
            if ischar(aTimeRange) && strcmp(aTimeRange, 'all')
                fFirstSample = 0;
                fLastSample  = Inf;
            else
                fFirstSample = aTimeRange(1) * aData.samplingFrequency;
                fLastSample  = aTimeRange(2) * aData.samplingFrequency;
            end
            
            fSparseChannelIndex = sparse(...
                double([aData.channelArray.channel.channelAchk] + 1) , ...
                double([aData.channelArray.channel.channelIndex] + 1), ...
                double(1:aData.channelArray.numChannels), ...
                BlockVectorLegacyLoader.MAX_CHANNEL_ARTICHOKE,...
                BlockVectorLegacyLoader.MAX_CHANNEL_INDEX);
            
            spikeIndex = 1;
            
            % Manage the spike array to avoid "growing" it
            fSpikeArrayCapacity = 0;
            
            fseek(aSourceData.FileID, int64(aSourceData.Start), 'bof');
            fStop = aSourceData.Start + aSourceData.EntryRecord.Length;
            
            while ftell(aSourceData.FileID) < fStop
                % Load the block header for this spike
                fCurrentSpike = [];
                fCurrentSpike.startingSample      = fread(aSourceData.FileID, 1, 'int64'); %Changed uint64 to int64
                
                if feof(aSourceData.FileID)
                    % we were at the end of the file
                    break;
                end
                
                fHardwareChannelIndex              = fread(aSourceData.FileID, 1, 'uint8');
                fHardwareChannelAchk               = fread(aSourceData.FileID, 1, 'uint8');
                fCurrentSpike.triggerSampleOffset  = fread(aSourceData.FileID, 1, 'uint32');
                fCurrentSpike.stDev                = fread(aSourceData.FileID, 1, 'double');
                fCurrentSpike.threshold            = fread(aSourceData.FileID, 1, 'double');
                
                % Get well and electrode row and column
                fCurrentSpike.channelArrayIndex = full(fSparseChannelIndex(fHardwareChannelAchk + 1, fHardwareChannelIndex + 1));
                if fCurrentSpike.channelArrayIndex == 0
                    error('load_AxIS_spike:invalidSpikeChannel', ...
                        ['Spike has invalid Artichoke/Channel number ' num2str(fHardwareChannelAchk) ...
                        ' / ' num2str(fHardwareChannelIndex)]);
                end
                
                fCurrentSpike.channelInfo = aData.channelArray.channel(fCurrentSpike.channelArrayIndex);
                
                % Historical channel number -- XY, where X = electrodeColumn and Y = electrodeRow
                % This field is deprecated, since it doesn't contain multiwell information.
                fCurrentSpike.channel = 10*fCurrentSpike.channelInfo.electrodeColumn + ...
                    fCurrentSpike.channelInfo.electrodeRow;
                
                % Seek forward to the end of the header.
                % It might seem cleaner to do this by saving the current position
                % before we read the header and seeking forward from there.
                % However, that method will not handle files larger than 4GB if fseek
                % and ftell are not 64-bit clean.
                fseek(aSourceData.FileID, aData.header.blockHeaderSize - Spike_v1.LOADED_HEADER_SIZE, 'cof');
                
                if ftell(aSourceData.FileID) >= fStop
                    % file ended before the start of the data -- looks like a truncated file
                    warning('load_AxIS_spike:headerTruncated', ...
                        'Spike header truncated at end of file');
                    break;
                end
                
                % Read the data for this spike
                fCurrentSpike.waveform = fread(aSourceData.FileID, aData.header.numSamplesPerBlock, 'int16=>int16');
                
                if length(fCurrentSpike.waveform) < aData.header.numSamplesPerBlock
                    % file ended in the middle of the data
                    warning('load_AxIS_spike:dataTruncated', ...
                        'Spike data truncated at end of file');
                    break;
                end
                
                % Does this spike match the channel filter?  If not, skip it.
                if ~ismember(fCurrentSpike.channelArrayIndex, aData.loadedChannels)
                    continue;
                end
                
                % Does this spike match the time filter?
                if fCurrentSpike.startingSample < fFirstSample
                    % Too early
                    continue;
                elseif fCurrentSpike.startingSample >= fLastSample
                    % Too late.  Since spikes are in chronological order, we
                    % can stop reading now.
                    break;
                end
                
                
                % Check for capacity in the spike array
                if spikeIndex > fSpikeArrayCapacity
                    if fSpikeArrayCapacity == 0
                        % Let it grow on its own
                        fSpikeArrayCapacity = 1;
                    else
                        % Double the capacity
                        fSpikeArrayCapacity = fSpikeArrayCapacity * 2;
                        aData.spikes(fSpikeArrayCapacity).startingSample = -1;
                    end
                end
                
                % save the spike
                aData.spikes(spikeIndex) = fCurrentSpike;
                spikeIndex = spikeIndex + 1;
            end
            
            % trim down the spike array
            if spikeIndex > 1
                aData.spikes = aData.spikes(1:spikeIndex-1);
            else
                aData.spikes = [];
            end
        end
    end
    
end