Fixes #141. Covers part of #147

+nf/extract.m

@@ -1,9 +1,9 @@
-%% Extract a specific subset of data from a neurofield output struct.
+%% Extract a specific subset of data from a nftsim output struct.
 %
 % The subset can be specified in terms of: traces; times; and nodes.
 %
 % ARGUMENTS:
-%        obj -- A neurofield output struct (a Matlab struct containing data
+%        obj -- A nftsim output struct (a Matlab struct containing data
 %               from a simulation).
 %        traces -- A cell array or a comma separated string of the traces
 %                  you want to extract, e.g. 'Propagator.2.phi, Coupling.2.nu'.
@@ -35,7 +35,7 @@
     % If no nodes are provided, output all nodes
     if ~isstruct(obj) || ~isfield(obj, 'data')
         error(['nf:' mfilename ':BadArgument'], ...
-              'The first argument must be a neurofield output struct.');
+              'The first argument must be a nftsim output struct.');
     end
 
     if nargin < 4 || isempty(nodes)

+nf/get_frequencies.m

@@ -2,20 +2,20 @@
 % Return corresponding grids of frequencies and positions
 %
 % ARGUMENTS:
-%        data -- .
-%        fs -- .
-%        Lx -- .
-%        Ly -- .
+%        data -- an 1D or 3D array of size space points along x, space points along y, tim points, 
+%        fs -- sampling temporal frequency in [Hz]
+%        Lx -- physical size of the spatial domain along x in [m] 
+%        Ly -- physical size of the spatial domain along y in [m]
 %
 % OUTPUT:
-%        f -- .
-%        Kx -- .
-%        Ky -- .
-%        x -- .
-%        y -- .
-%        df -- .
-%        dk -- .
-%        dx -- .
+%        f -- vector of (positives) frequencies in [Hz]
+%        Kx -- vector of angular wavenumbers in [rad/m]
+%        Ky -- vector of angular wavenumbers in [rad/m]
+%        x -- vector with spatial coordinates in [m]
+%        y -- vector with spatial coordinates in [m]
+%        df -- temporal frequency resolution in [Hz]
+%        dk -- radial spatial frequency resolution in [rad/m]
+%        dx -- spatial resolution in [m]
 %
 % REQUIRES:
 %           -- <description>
@@ -38,16 +38,19 @@
     x = [];
     y = [];
 
+    % This is the single positive sided frequency vector
     if isvector(data)
-        f = (0:(fs / length(data)):(fs / 2)).'; % This is the single sided frequency
+        f = (0:(fs / length(data)):(fs / 2)).'; 
         return
     end
 
-    f = (0:(fs / size(data, 3)):(fs / 2)).'; % This is the single sided frequency
+    % This is the single positive sided frequency
+    f = (0:(fs / size(data, 3)):(fs / 2)).'; 
     df = fs / size(data, 3);
-
-    Kx = 2*pi*(0:1/Lx:size(data, 1)/Lx/2); % This is the single sided frequency
-    Ky = 2*pi*(0:1/Ly:size(data, 2)/Ly/2); % This is the single sided frequency
+
+    % Wavenumbers / single sided and positive
+    Kx = 2*pi*(0:1/Lx:size(data, 1)/Lx/2);
+    Ky = 2*pi*(0:1/Ly:size(data, 2)/Ly/2);
 
     if mod(size(data, 1), 2) % If there is NO nyquist frequency component
         Kx = [-Kx(end:-1:2) Kx];
@@ -60,13 +63,24 @@
     else
         Ky = [-Ky(end:-1:2) Ky(1:end-1)];
     end
+
     [Kx, Ky] = meshgrid(Kx, Ky);
+
+    % Smallest wavenumber - resolution in k-space (angular spatial frequency)
+    % This is assuming the spatial domain is square (Lx=Ly)
+
     dk = 2 * pi / Lx;
+
+    % dx and dy should be the same
+    dx = Lx / size(data, 1); % [m]
+    dy = Ly / size(data, 1); % [m]
+
+    % the stencil in NFTsim assumes the coordinates of a parcel of the discretized domain
+    % is at the centre of the parcel, thus the first coordinate is at (dx/2; dy/2) 
+
+    x  = dx * (0:size(data, 1)) + dx/2;
+    y  = dy * (0:size(data, 1)) + dx/2;
 
-    dx = Lx / size(data, 1); % Metres per pixel
-    x = dx * (0:size(data, 1)-1);
-    dy = Ly / size(data, 1);
-    y = dy * (0:size(data, 1)-1);
     [x, y] = meshgrid(x, y);
 
 end %function get_frequencies()

+nf/grid.m

@@ -8,7 +8,7 @@
 % output all nodes), or the number of nodes is not a perfect square..
 %
 % ARGUMENTS:
-%        obj   -- A neurofield output struct (a Matlab struct containing data
+%        obj   -- A nftsim output struct (a Matlab struct containing data
 %                 from a simulation).
 %        trace -- A string with the name of the array to reshape. 
 %
@@ -35,7 +35,7 @@
 
     if output_nodes ~= obj.input_nodes 
         error(['nf:' mfilename ':IncompatibleOutput'], ...
-              'Output from NeuroField must be for all nodes')
+              'Output from NFTsim must be for all nodes')
     end    
 
     data = nf.extract(obj, trace);
@@ -48,7 +48,7 @@
         longside_nodes  = obj.longside_nodes;
         shortside_nodes = obj.input_nodes / obj.longside_nodes;
     end
-    %Reshape to an array of (n,m,tpts)
+    %Reshape to an array of (nx,ny,tpts)
 
     %data = reshape(data, grid_edge, grid_edge, obj.npoints);
 

+nf/plot_timeseries.m

@@ -17,7 +17,7 @@
 %    figure_handles -- cell array of figure handles.
 %
 % REQUIRES:
-%    nf.extract() -- Extract a specific subset of data from a neurofield
+%    nf.extract() -- Extract a specific subset of data from a nftsim
 %                    output struct.
 %
 % AUTHOR:
@@ -28,7 +28,7 @@
 %{
     %Either run a simulation:
     nf_obj = nf.run('./configs/eirs-corticothalamic.conf')
-    %Or load some neurofield output data
+    %Or load some nftsim output data
     nf_obj = nf.read('./configs/eirs-corticothalamic.output')
 
     %Plot every fourth node for the trace 'Propagator.1.phi'.

+nf/read.m

@@ -1,11 +1,11 @@
-%% Read a neurofield output file and return a neurofield output struct.
+%% Read a nftsim output file and return a nftsim output struct.
 %
 % ARGUMENTS:
-%        fname -- The name of the neurofield output file to read (absolute
+%        fname -- The name of the nftsim output file to read (absolute
 %                 or relative path).
 %
 % OUTPUT:
-%        obj -- A neurofield output struct. A Matlab struct containing data
+%        obj -- A nftsim output struct. A Matlab struct containing data
 %               and parameters from a simulation.
 %
 % AUTHOR:
@@ -42,7 +42,7 @@
 
 
     while isempty(strfind(buffer, '======================='))
-        % TODO: consider cleaning up this part.
+        % TODO: CLEAN UP - this part refers to EEGCODE.
         if ~isempty(strfind(buffer, 'Time  |'))
             error(['nf:' mfilename ':OldStyleOutput'], ...
                   'Did you try and open and old-style output file? Found a | that looked like a delimiter.')

+nf/report.m

@@ -1,7 +1,7 @@
-%% Given a neurofield output struct, print some information about it.
+%% Given a nftsim output struct, print some information about it.
 %
 % ARGUMENTS:
-%        obj -- A neurofield output struct (a Matlab struct containing data
+%        obj -- A nftsim output struct (a Matlab struct containing data
 %               from a simulation).
 %
 % OUTPUT: Prints to terminal.
@@ -16,7 +16,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 function report(obj)
-    fprintf(1, 'Output for: neurofield -i "%s"\n', obj.conf_file)
+    fprintf(1, 'Output for: nftsim -i "%s"\n', obj.conf_file)
     fprintf(1, 'Traces: ');
     for j = 1:length(obj.fields)
         fprintf(1, '%s ', obj.fields{j});

+nf/rfft.m

@@ -21,6 +21,7 @@
 %        http://www.brain-dynamics.net/~chris_rennie/fourier.pdf
 %
 % AUTHOR:
+%     Original: Chris Rennie circa 2000
 %     Romesh Abeysuriya (2012-03-22).
 %
 % USAGE:
@@ -75,9 +76,10 @@
         P(2:(end - 1), :) = P(2:(end - 1), :) .* 2.0;
     end
 
+    % TODO: check normalizations
     % P now contains properly normalized spectral power
     %f = fs / 2 * linspace(0, 1, NFFT / 2 + 1)';
-    f = 0:(fs / NFFT):(fs / 2); % I think this is more correct
+    f = 0:(fs / NFFT):(fs / 2); % I think this is more correct 
     P = mean(P, 2);
     P = P ./ f(2); % Divide by frequency bin size to get power density
 

+nf/run.m

@@ -1,38 +1,38 @@
-%% Function to run neurofield and return a neurofield output struct.
+%% Function to run nftsim and return a nftsim output struct.
 %
-% Provided a configuration file-name (fname.conf), run the neurofield
+% Provided a configuration file-name (fname.conf), run the nftsim
 % executable, generating an output file (fname.output). Optionally, if an
 % output argument is provided then, parse the output file and return a
-% neurofield output struct containing the simulation results.
+% nftsim output struct containing the simulation results.
 %
 %
 % ARGUMENTS:
 %        fname -- Name of the configuration file, it can be with or without
 %                 the .conf extension.
 %        time_stamp -- boolean flag to use a time_stamp YYYY-MM-DDTHHMMSS
 %                      in the output file name.
-%        neurofield_path -- neurofield executable (full or relative path).
+%        nftsim_path -- nftsim executable (full or relative path).
 %
 % OUTPUT: Writes a .output file in the same location as the .conf file.
-%        obj -- A neurofield output struct (a Matlab struct containing
+%        obj -- A nftsim output struct (a Matlab struct containing
 %               data from a simulation).
 %
 % REQUIRES:
-%        neurofield -- The neurofield executable, must be in your path.
-%        nf.read -- Read a neurofield output file and return a neurofield
+%        nftsim -- The nftsim executable, must be in your path.
+%        nf.read -- Read a nftsim output file and return a nftsim
 %                   output struct.
 %
 % AUTHOR:
 %     Romesh Abeysuriya (2012-03-22).
 %
 % USAGE:
 %{
-    %At a matlab command promt, from neurofield's main directory:
+    %At a matlab command promt, from nftsim's main directory:
     nf_obj = nf.run('./configs/eirs-corticothalamic.conf')
 %}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function obj = run(fname, time_stamp, neurofield_path)
+function obj = run(fname, time_stamp, nftsim_path)
     %
     tic;
     fname = strrep(fname, '.conf', ''); %Strip any .conf suffix.
@@ -41,39 +41,39 @@
     if nargin < 2
         time_stamp = false;
     end
-    % If we were not provided a path to neurofield, try to determine one.
-    if nargin < 3 || isempty(neurofield_path)
+    % If we were not provided a path to nftsim, try to determine one.
+    if nargin < 3 || isempty(nftsim_path)
         % Check typical locations, the first path that exists will be selected.
-        locations = {'neurofield', ...
-                     './bin/neurofield', ...
-                     './neurofield/bin/neurofield', ...
-                     'neurofield.exe'};
+        locations = {'nftsim', ...
+                     './bin/nftsim', ...
+                     './nftsim/bin/nftsim', ...
+                     'nftsim.exe'};
         selected_path = find(cellfun(@(name) exist(name, 'file')==2, locations), 1, 'first');
         if ~isempty(selected_path)
-            neurofield_path = locations{selected_path};
+            nftsim_path = locations{selected_path};
         else
             error(['nf:' mfilename ':BadPath'], ...
-                  'neurofield not found. Either change into the neurofield base directory or make a symlink to neurofield in the current directory.');
+                  'nftsim not found. Either change into the nftsim base directory or make a symlink to nftsim in the current directory.');
         end
     % If we were provided a path, check that it is valid.
-    elseif ~exist(neurofield_path, 'file')
+    elseif ~exist(nftsim_path, 'file')
         error(['nf:' mfilename ':BadPath'], ...
-              'The neurofield_path you provided is incorrect:"%s".',neurofield_path);
+              'The nftsim_path you provided is incorrect:"%s".',nftsim_path);
     end
 
     if ~time_stamp
-        neurofield_cmd = sprintf('%s -i %s.conf -o %s.output', neurofield_path, fname, fname);
+        nftsim_cmd = sprintf('%s -i %s.conf -o %s.output', nftsim_path, fname, fname);
     else
-        neurofield_cmd = sprintf('%s -i %s.conf -t', neurofield_path, fname);
+        nftsim_cmd = sprintf('%s -i %s.conf -t', nftsim_path, fname);
     end
     fprintf('INFO: Executing command:\n')    
-    fprintf('%s\n', neurofield_cmd);
-    [status, cmdout] = system(neurofield_cmd);
+    fprintf('%s\n', nftsim_cmd);
+    [status, cmdout] = system(nftsim_cmd);
 
     string(cmdout)
     if status ~= 0
-        error(['nf:' mfilename ':NeurofieldError'], ...
-              'An error occurred while running neurofield!');
+        error(['nf:' mfilename ':NFTsimError'], ...
+              'An error occurred while running nftsim!');
     end
 
     fprintf(1, 'INFO: tic-toc: Took about %.3f seconds\n', toc);

+nf/spatial_spectrum.m

@@ -10,19 +10,19 @@
 %        spatial_filter -- set to 1 to enable the usual exponential k filter
 %
 % OUTPUT:
-%        f -- .
-%        P -- .
-%        Kx -- .
-%        Ky -- .
-%        Pkf -- .
-%        x -- .
-%        y -- .
-%        Prf -- .
+%        f -- vector of temporal frequencies in [Hz]
+%        P -- power spectral density in [au / Hz ???]????
+%        Kx -- vector with angular wavenumbers in [rad/m]
+%        Ky -- vector with angular wavenumbers in [rad/m]
+%        Pkf -- power spectral density? array in [units?????]
+%        x -- array of spatial coordinates along x [m]
+%        y -- array of spatial coordinates along y [m]
+%        Prf --power spectral density ? in space frequency [units?????]
 %
 % REQUIRES:
-%        nf.get_frequencies() -- <description>
-%        nf.grid()  -- <description>
-%        nf.partition() -- <description>
+%        nf.get_frequencies() 
+%        nf.grid()  
+%        nf.partition() 
 %
 % REFERENCES:
 %
@@ -119,9 +119,9 @@
 
 
 function [Pkf, Prf] = get_spectrum(data, k_mask, k_filter, Lx, fs)
-    df = fs / (size(data, 3));
-    dk = 2 * pi / Lx;
-    dx = Lx / size(data, 1); % Metres per pixel
+    df = fs / (size(data, 3)); % [Hz]
+    dk = 2 * pi / Lx;          % spatial frequency resolution in [rad/m] -- smallest nonzero wavenumber
+    dx = Lx / size(data, 1);   % [m]
 
     % First, get the 3D FFT and apply volume conduction
     P = fftshift(fftn(data));
@@ -155,7 +155,7 @@
             Prf(:, :, 2:(end - 1)) = Prf(:, :, 2:(end - 1)) * 2;
         end
 
-        % Divide by N, dx and dx so that integration gives the correct power
+        % Divide by N, dx and dx so that integration gives the correct power or PSD?????
         % Division by N on the grounds that energy in real space is given by summing abs(x)^2/N
         % So we could call the power P=abs(x)^2/N so that sum(P) gives the correct power
         % And then division by dx*dx means that integral(P) over position gives the right answer