diff --git a/aod/AodHRVolumeReader.m b/aod/AodHRVolumeReader.m
new file mode 100755
index 0000000..5f284ec
--- /dev/null
+++ b/aod/AodHRVolumeReader.m
@@ -0,0 +1,114 @@
+
+classdef AodHRVolumeReader < HDF5Helper
+    properties
+        x = [];
+        y = [];
+        z = [];
+        reps = [];
+        version = [];
+        vfilename = [];
+    end
+    
+    methods
+        function self = AodHRVolumeReader(filename,vfilename)
+            
+            self = self@HDF5Helper(filename, 'ImData');
+            
+            isHRV = (H5Tools.existAttribute(self.fp,'CellFocussedVolume') && ...
+                H5Tools.readAttribute(self.fp, 'CellFocussedVolume') == 1);
+            self.version = H5Tools.readAttribute(self.fp, 'Version');
+            assert((H5Tools.existAttribute(self.fp,'ScanType') && ...
+                strncmp(H5Tools.readAttribute(self.fp, 'ScanType')', 'Volume', 6) == 1) ...
+                || isHRV , ...
+                'Not a volume scan file');
+            
+            
+            coord = AodCoordinatesReader(self.filename,self.version);
+            coords =  coord(:);
+            numPoints = size(coords,1);
+            self.reps = floor(self.sz(1) / numPoints);
+            
+            ar = aodReader(vfilename,'Volume');
+            
+            self.x = ar.x;
+            self.y = ar.y;
+            self.z = ar.z;
+            self.vfilename = vfilename;
+            
+        end
+        
+        function volume = subsref(self, s)
+            
+            ar = aodReader(self.filename,'Functional');
+            data = ar(:,:,:);
+            coord = AodCoordinatesReader(self.filename,self.version);
+            coords =  coord(:);
+            volume = nan(length(self.x),length(self.y),length(self.z),size(data,1),size(data,3));
+            
+            for i = 1:size(coords,1)
+                xind = find(single(coords(i,1))==single(self.x));
+                yind = find(single(coords(i,2))==single(self.y));
+                zind = find(single(coords(i,3))==single(self.z));
+                volume(xind,yind,zind,:,:) = squeeze(data(:,i,:));
+            end
+            
+            if(strcmp(s(1).type,'()') == 0)
+                if 0 && ismember(s.subs,properties(self))
+                    volume = self.(s.subs);
+                    return;
+                end
+                volume = subsref@HDF5Helper(self,s);
+                return
+            end
+            
+            % make sure subscripting has the right form
+            assert(numel(s) == 1 && strcmp(s.type, '()') && numel(s.subs) == 5 || ...
+                numel(s) == 1 && strcmp(s.type, '()') && numel(s.subs) == 4 || ...
+                numel(s) == 1 && strcmp(s.type, '()') && numel(s.subs) == 3, ...
+                'MATLAB:badsubscript', 'Only subscripting of the form (samples, channels) is allowed!')
+            
+            % samples and channels
+            x = s(1).subs{1};
+            y = s(1).subs{2};
+            z = s(1).subs{3};
+            
+            if numel(s(1).subs)<5
+                chan = 1;
+            else
+                chan = s(1).subs{5};
+            end
+            
+            if numel(s(1).subs) < 4
+                reps = 1:self.reps;
+            else
+                reps = s(1).subs{4};
+            end
+            
+            % all samples requested?
+            if ~HDF5Helper.iscolon(x)
+                volume = volume(x,:,:,:,:);
+            end
+            
+            % all samples requested?
+            if ~HDF5Helper.iscolon(y)
+                volume = volume(:,y,:,:,:);
+            end
+            
+            % all samples requested?
+            if ~HDF5Helper.iscolon(z)
+                volume = volume(:,:,z,:,:);
+            end
+            
+            % all samples requested?
+            if ~HDF5Helper.iscolon(reps)
+                volume = volume(:,:,:,reps,:);
+            end
+            
+            % all channels requested?
+            if ~HDF5Helper.iscolon(chan)
+                volume = volume(:,:,:,:,chan);
+            end
+            
+        end
+    end
+end
\ No newline at end of file
diff --git a/aod/AodScanReader.asv b/aod/AodScanReader.asv
deleted file mode 100755
index 97f6874..0000000
--- a/aod/AodScanReader.asv
+++ /dev/null
@@ -1,133 +0,0 @@
-
-classdef AodScanReader < HDF5Helper
-    properties
-        coordinates = [];
-        numPoints = 0;
-        reshapedSize = [];
-        motionIncrement = 0;
-    end
-    
-    methods 
-        function self = AodScanReader(filename)
-            %self.restrict()
-
-            self = self@HDF5Helper(filename, 'ImData');
-           
-            version = H5Tools.readAttribute(self.fp,'Version');
-            
-            %%%% Coordinates are now stored as data 2014-12-10 %%%%%%%%%%%%
-            try
-                coord = AodCoordinatesReader(filename,version);
-                self.coordinates = coord(:);
-            catch
-                self.coordinates = reshape(H5Tools.readAttribute(self.fp, ...
-                    'PointCoordinates'),3,[])';
-                self.coordinates = bsxfun(@times, self.coordinates, ...
-                  [1/1460000 1/1460000 1/700]); % [acq.AodScan.x_step acq.AodScan.y_step acq.AodScan.z_step]);
-            end
-          
-            
-            
-            self.numPoints = size(self.coordinates,1);
-            timepoints = floor(self.sz(1) / self.numPoints);
-            
-            if H5Tools.existAttribute(self.fp, 'MotionIncrement')
-                self.motionIncrement = H5Tools.readAttribute(self.fp ,'MotionIncrement');
-                self.Fs = self.Fs / self.numPoints  * (self.motionIncrement - 1) / self.motionIncrement;
-            else
-                self.Fs = self.Fs / self.numPoints;
-            end
-            
-            self.reshapedSize = [timepoints self.numPoints self.sz(2)];
-        end
-        
-        function [x, vargout] = size(self, s)
-            x = self.reshapedSize;
-            if nargin > 1
-                x = x(s);
-            end
-            
-            if nargout == 1
-                % unchanged
-            elseif nargout == 2
-                vargout = x(2);
-                x = x(1);
-            end
-        end
-        
-        function x = subsref(self, s)
-            % make sure subscripting has the right form
-            
-            if(strcmp(s.type,'()') == 0) 
-                if 0 && ismember(s.subs,properties(self))
-                    x = self.(s.subs);
-                    return;
-                end
-                x = subsref@HDF5Helper(self,s);
-                return
-            end
-            
-            assert(numel(s) == 1 && strcmp(s.type, '()') && numel(s.subs) == 3 || ...
-                numel(s) == 1 && strcmp(s.type, '()') && numel(s.subs) == 2, ...
-                'MATLAB:badsubscript', 'Only subscripting of the form (samples, channels) is allowed!')
-            
-            % samples and channels
-            samples = s(1).subs{1};
-            points = s(1).subs{2};
-
-            % all samples requested?
-            if HDF5Helper.iscolon(samples)
-                samples = 1:self.reshapedSize(1);
-            else
-                % Check for valid range of samples
-                assert(all(samples <= self.reshapedSize(1) & samples > 0), 'MATLAB:badsubscript', ...
-                    'Sample index out of range [1 %d]', self.reshapedSize(1));
-            end
-
-            % time channel requested?  Shortcut out
-            if ischar(points) && points == 't'
-                assert(self.t0 > 0, 't0 has not been updated in this file!')
-                if HDF5Helper.iscolon(samples)
-                    x = self.t0 + 1000 * (0:self.reshapedSize(2)-1)' / self.Fs;
-                else
-                    x = self.t0 + 1000 * (samples(:)-1)' / self.Fs;
-                end
-                return;
-            end
-
-            if numel(s(1).subs) >= 3
-                channels = s(1).subs{3};
-            else 
-                channels = 1;
-            end
-            
-            % all points requested?
-            if HDF5Helper.iscolon(points)
-                points = 1:self.reshapedSize(2);
-            else
-                % Check for valid range of channels
-                assert(all(points <= self.reshapedSize(2) & points > 0), ...
-                    'MATLAB:badsubscript', 'Point index out of range [1 %d]', self.reshapedSize(2));
-            end
-            
-            % all channels requested?
-            if HDF5Helper.iscolon(channels)
-                channels = 1:self.reshapedSize(3);
-            else
-                % Check for valid range of channels
-                assert(all(channels <= self.reshapedSize(3) & channels > 0), ...
-                    'MATLAB:badsubscript', 'Channel index out of range [1 %d]', self.reshapedSize(3));
-            end
-
-            i1 = repmat(samples,numel(points),1);
-            i2 = repmat(points,1,numel(samples));
-            
-            ind = sub2ind(self.reshapedSize([2 1]), i2(:), i1(:));
-            
-            s.subs = {ind, channels};
-            x = subsref@HDF5Helper(self,s);
-            x = reshape(x,[numel(points) numel(samples) numel(channels)]);
-            x = permute(x,[2 1 3]);
-        end
-    end
-end
\ No newline at end of file
diff --git a/aod/aodReader.m b/aod/aodReader.m
index 1b69e58..d17c5e3 100755
--- a/aod/aodReader.m
+++ b/aod/aodReader.m
@@ -45,6 +45,8 @@
     br = AodMotionReader(fileName);
 elseif strcmp(dataset, 'Volume') == 1
     br = AodVolumeReader(fileName);
+elseif strcmp(dataset, 'HRVolume') ==1
+    br = AodHRVolumeReader(fileName,varargin{1});
 else
     error('Unknown dataset.  Options are Temporal, Functional,Volume and Motion.');
 end
diff --git a/raw/HDF5Helper.m b/raw/HDF5Helper.m
new file mode 100755
index 0000000..697dadb
--- /dev/null
+++ b/raw/HDF5Helper.m
@@ -0,0 +1,150 @@
+classdef HDF5Helper
+    properties
+        filename = '';
+        fp = 0;
+        sz = [];
+        dataset = '';
+        t0 = 0;
+        Fs = 0;
+    end
+    
+    methods
+        function self = HDF5Helper(filename, dataset)
+            self.filename = filename;
+            self.fp = H5Tools.openFamily(self.filename);
+            self.dataset = dataset;
+            self.sz = fliplr(H5Tools.getDatasetDim(self.fp, self.dataset));
+
+            if(H5Tools.existAttribute(self.fp, 't0'))
+                self.t0 = H5Tools.readAttribute(self.fp, 't0');
+            else
+                self.t0 = 0;
+            end
+
+             if(H5Tools.existAttribute(self.fp, 'Fs'))
+                self.Fs = H5Tools.readAttribute(self.fp, 'Fs');
+             else
+                self.Fs = 50000;
+             end
+        end
+        
+        function self = close(self)
+            if ~isempty(self.fp)
+                H5F.close(self.fp);
+                self.fp = [];
+            end
+        end
+
+        function Fs = getSamplingRate(self)
+            % Return the sampling rate
+            Fs = self.Fs;
+        end
+        
+        function l = length(self)
+            l = self.sz(1);
+        end
+
+        function x = end(self,dim,varargin)
+            x = size(self,dim);
+        end
+
+        function varargout = size(self, dim)    
+            if (nargout == 1) && (nargin > 1)
+                varargout{1} = self.sz(dim);
+            elseif (nargout > 1)
+                varargout = num2cell( self.sz(1:nargout) );
+            else
+                varargout{1} = self.sz;
+            end
+        end
+        
+        function x = subsref(self, s)
+            % Subscripting.
+            %   x = br(samples, channels). channels can be either channel indices or
+            %   't' for the timestamps in milliseconds.
+            %
+            % AE 2011-04-11
+            
+            if (strcmp(s(1).type, '.') == 1)
+                assert(ismember(s(1).subs,properties(self)), 'MATLAB:badsubscript', 'Field not found');
+                x = self.(s(1).subs);
+                
+                if length(s) > 1
+                    x = subsref(x,s(2:end));
+                end
+                
+            	return;
+            end
+            
+            % make sure subscripting has the right form
+            assert(numel(s) == 1 && strcmp(s.type, '()') && numel(s.subs) == 2, ...
+                'MATLAB:badsubscript', 'Only subscripting of the form (samples, channels) is allowed!')
+            
+            % samples and channels
+            samples = s(1).subs{1};
+            channels = s(1).subs{2};
+            
+            % all samples requested?
+            if HDF5Helper.iscolon(samples)
+                nSamples = self.sz(1);
+            else
+                % Check for valid range of samples
+                assert(all(samples <= self.sz(1) & samples > 0), 'MATLAB:badsubscript', ...
+                    'Sample index out of range [1 %d]', self.sz(1));
+                nSamples = numel(samples);
+            end
+            
+            % time channel requested?
+            if ischar(channels) && channels == 't'
+                assert(self.t0 > 0, 't0 has not been updated in this file!')
+                if iscolon(samples)
+                    x = self.t0 + 1000 * (0:self.sz(1)-1)' / self.Fs;
+                else
+                    x = self.t0 + 1000 * (samples(:)-1)' / self.Fs;
+                end
+            else
+                
+                % all channels requested?
+                if HDF5Helper.iscolon(channels)
+                    channels = 1:(self.nbImChannels);
+                else
+                    % Check for valid range of channels
+                    assert(all(channels <= self.sz(2) & channels > 0), ...
+                        'MATLAB:badsubscript', 'Channel index out of range [1 %d]', self.sz(2));
+                end
+                nChannels = numel(channels);
+                
+                % Convert to actual channel numbers in the recording file
+                %channels = self.imChIndices(channels);
+                
+                x = zeros(nSamples, nChannels);
+                
+                if HDF5Helper.iscolon(samples)
+                    % reading all samples
+                    for i = 1:nChannels
+                        x(:,i) = H5Tools.readDataset(self.fp, self.dataset, 'range', [channels(i), 1], [channels(i), self.sz(1)]);
+                    end
+                elseif length(samples) > 2 && samples(end) - samples(1) == length(samples) - 1 && all(diff(samples) == 1)
+                    % reading continuous block of samples
+                    for i = 1:nChannels
+                        x(:,i) = H5Tools.readDataset(self.fp, self.dataset, 'range', [channels(i), samples(1)], [channels(i), samples(end)]);
+                    end
+                else
+                    % reading arbitrary set of samples
+                    for i = 1:nChannels
+                        x(:,i) = H5Tools.readDataset(self.fp, self.dataset, 'index', [repmat(channels(i), nSamples, 1) , samples(:)]);
+                    end
+                end
+                
+                % scale to (micro/milli?)volts
+            end
+        end
+    end
+    
+    methods(Static)
+        function b = iscolon(x)
+            b = ischar(x) && isscalar(x) && x == ':';
+        end
+    end
+
+end
\ No newline at end of file
diff --git a/setPath.m b/setPath.m
index 6e3ef5c..adb5f7c 100755
--- a/setPath.m
+++ b/setPath.m
@@ -1,6 +1,7 @@
 function setPath
     % Add subfolders for HDF5 access to the path
     addpath(fileparts(mfilename('fullpath')))
+    addpath(fullfile(fileparts(mfilename('fullpath')), 'aod'))
     addpath(fullfile(fileparts(mfilename('fullpath')), 'raw'))
     addpath(fullfile(fileparts(mfilename('fullpath')), 'spikes'))
 end