MATLAB FUNCTIONS FOR MOTOR PSYCHOPHYSICS:
(c) Authors:     Moshe Abeles (M.A.)  <abeles@md.huji.ac.il>
                             Vladimir Litvak (V.L.) <vladimir@md2.huji.ac.il>
                             Yoram Ben Shaul (Y.B.S.) <ybss@md2.huji.ac.il>
                             Zoltan Nadasdy (Z.N.) <zoli@md2.huji.ac.il>
Keys:
bold: function name
red: obsolete

(..) Parameter field.
[..] key stroke


function all2ascii()
to convert movement object to ascii list
of X,Y

was made for shipping data to TAMAR FLASH



function all2new(dirname)


function analyse(fileName)


script animatemove
DOES: ANIMATEMOVE a script to animate a movement.
this will show first the motion, then the motion +tangvelocity,
then the motion + tangvelocity + powerlaw fit.
this is a script that expects the following parameters to
be ready:
x
y
t
tv
fileName
theta
rewardTimes
segFile
movFile
resFile
see PREPAREANIMATE for how to extact these data
  set DELAY to suit your machine

Jan-1999  MA
Updates
   May-1999  set the SHOW_TRGTS flag if you want to see the targets
             together with the animation.
             Q after inquiery will quit the program



function automatic()


function automatic_inter_reward()
DOES:    The function parses all the 'digital' files
found in the curses all the 'digital' files
found in the current directory, determines which
of them contain scribling trials, calculates inter reward
intervals in these files, sorts the results by dates and
plots median inter reward interval vs. date.


function box3d (x,y,Z,MC)
DOES:  box3d  to display a 3d histogram in box fashion
USE: [f3d, ax3d] = box3d (x,y,Z,MC)
INPUT: x and y contain the baseline coordinates
Z contains the count in each bin
MC is the colormap arrey to be used
returns a handle to the 3d figure in f3d, and
to the axes in ax3d.

temporary generate the figure
UPDATES:
  Feb-1996 fitted for show_loc_freq
  16-Aug-1996  colorbar is NOT plotted  M.A.
  17-Mar-1997  ColorMap added, rotation angle improved

function [f3d, ax3d] = box3d (x,y,Z,MC)

find the largest existing figure



script buzztry

Beeps when the hand trajectory hits the target.



script delfig
DOES: DELFIG to find an object with some property
and if exists delete it.
e.g. delFig('Tag','TangVel')

Jan 1999 M.A.


function experiment2mat(inFile, outDir, digFile)
EXPERIMENT2MAT convert analog file to movement object
and save it to a matlab file.
EXPERIMENT2MAT(INFILE, OUTDIR, DIGDIR)
INPUT: INFILE the analog file na
INPUT: INFILE the analog file name (without a directory)
OUTDIR where to put the MAT file
DIGFILE name of event file
this procedure uses 2 scritps: MOV2MAT and GETXY
UPDATES:
Jan-1999
  times of rewards is read from event file and are added
  to movement in cell arrey M.EVENTS
  if outDir exists the MAT file is put there
  The time of first sample in the file is added.  MA
  Either <xh,yh> or <HandX,HandY> are used according to
   the one with smaller jitter in the tangential velocity.  MA
Feb 1999
Added to the cell array m.events the qualifier of the relevant trial.
This change is dependent on changes made in parse_events_file(). YBS

Since encoder was giving wrong readings- i use only the readings of the
BC - YBS.



function fig2eps (fileName)
DOES: to convert a figure to a COLOR eps fert a figure to a COLOR eps file with
TIFF preview.
example: click on the figure you wish to save
  then fig2eps('pictures\fig1.eps')

Mar-1999 MA



function filt_mov (hx, dt, f0)
USE: function shx  = filt_mov (hx, dt, f0)
FILT_MOV Filters input data by low-pass.
FILT_MOV(data, dt, cutoff)
DOES: Filtering is done by Gaussian transfer function (exp(-(f/cutoff)^2/2)
in the freq. domain
hx    - original data
cutoff - the std-dev in the freq. domain
shx = the smoothed data.
WARNING:  The first few points and last few points are
grossly displaced!
Superficial testing with scribling shows that f0=4Hz
is great!

Feb 1998  MA
UPDATES: Modified 1/9/98 by Vladimir Litvak for use in movement toolbox.
Modified 11/10/98 by Vladimir Litvak - mirror image of movement added to
    prevent distortion of the movement by fft.


function findpeaks(data, filt, thresh, absthresh)
USE: function [startMinPoints, endMinPoints , maxPoints]=findpeaks(data, filt, thresh, absthresh)
DOES: FINDPEAKS Finds significant minimal and maximal points in the data.
The maximal points are chosen to be higher than a certain threshold
thresh, which is defined in terms of the mean of the data. For instance,
if thresh==3/4 peaks are defined as local maxima higherre defined as local maxima higher than
3/4*mean(data). Absolute thhreshold value can also be defined by using
findpeaks(data, filt, thresh, 'abs').
The minimal points are low local minima which separate two adjacent peaks.
To avoid finding too many local minima the data is smoothed at the first
stage. filt parameter defines the degree of smoothing - filt neighboring
points are averaged.
Default values are thresh==2/3, filt==8.

AUTHOR: Written by Vladimir Litvak Aug-98
UPDATES: Modified 1/9/98 for use in movement toolbox.
WARNING
filt must be even!

UPDATES:
  Jan-1999
    filt is forced to be even
    default values for <filt,thresh> were changed from <12,1/2>
        to <4,1/3>.   MA


function fiterror(mov, params)
USE: function err=fiterror(mov, params)
DOES: FITERROR Calculates average least squares error of pwrlawfit
err=PWRLAWFIT(mov, params)
INPUT: mov - movement under study.
params - parameters of the fit.

AUTHOR: Written by Vladimir Litvak 20/10/98


function genpwrlaw(params, xdata)
function v=genpwrlaw(params, xdata)

function get_s_e(Si,Ei,absVh)
USE:  Get arreys for pre-rest, motion, post rest.
  index of 0 - means does not exist
&nbsf 0 - means does not exist
  call by [eb , mv , ea] = get_s_e(Si,Ei,absVh)

INPUT:
  absVh  arrey with absolute velocities
  Si    , Ei    arreys of indices where motion starts/ends

OUTPUT:
  eb  Indices of rest before a motion starts
  mv  Indices of the motion itself
  ea  Indices of rest after the end of motion
  size of all three output matrices is (m,2) where
    m is the number of movements, the first column is the
    start index and the second is the end index.

function [eb , mv , ea] = get_s_e(Si,Ei,absVh)



function get_files(fileNumbers)


script getManip
DOES:   A script to compute the manipulandum data call by  getManip;
  calling should be done only after  fname  was assigned
   e.g:
      >> fname = '03aa_050';
      >> getManip

   where fname is a name of an analog data file as generated
   by ana2mat, without the .mat extention.
   getManip will generate/modify the following:
   HandX , HandY - the (x,y) coordinates in mm as read from the file
   alpha , beta  - the (upperalpha , beta  - the (upper,lower) arm encoders in degrees, as read
                   from the file.
   X_accel , Y_accel - the accelerometer reading (unknown units) , as
                       read from the file, after subtracting the amplifier
                       offsets at zero (as found in files with no motion).
   X_accel1 , Y_accel1 - temporary accelerations after scaling to fit to
                         units of mm/sec^2 (see getAccel).
   HandXdt , HandYdt , ... the sampling interval in seconds.
   b_rad - the angle of the lower arm relative to the X-axiss, in radians.
   xh , yh - the (x,y) cordinates as computed from (alpha,beta).
   Addx , Addy - the accelerations in hand coordinates, as
                 reconstructed from (X_accel,Y_accel).
   Adx , r>   Adx , Ady   - velocities in Kand coordinates as computed from Addx,Addy
   Ax_coef, Ay_coef - temporary for fitting accelerations.
   Ei - list of indices where motion stopped (see whenMov.m)
   Hax   Hay   - accelerations as computed from (xh,yh) and converted
                 to accelerometer coordinates.
   Hddx , Hddy - accelerations as computed from (xh,yh).
   Hdx   Hdy   - velocities as computed from (xh,yh).
   Si - list of indices where motion started (see whenMov.m)
   X_coef  Y_coef - temporary, coeficients for scaling X_accel (Y_accel)
                    based on fitting of Hax to X_accel (Hay to Y_accel).
   Ax_coef , Ay_coef - temporary, coeficients for scaling X_accel (Y_accel)
                    based on fitting of X_accel to Hddx(Y_accel to Hddy).
   absAh - the absolute acceleration (based on xh,yh) in mmatlab-mode.elm/sec^2.
   absAcc - the absoli>absAcc
- the absolute acceleration computed from (X_accel,Y_accel)
   absVh - the absolute velocity (in mm.sec) computed from (xh,yh).
   moving - 0 when NOT moving 1 when IS moving.
   ns - no of samples in this file
   run - minimal no of 0's in a row (in MOVING) to be consideres
         as a real stop of motion
   th  - threshold for moving
Summary of conversions

                         |Handle coord.   | Accel. coord.      |
-------------------------|----------------|--------------------|
from (xh , yh)           |Hddx, Hddy      | Hax, Hay           |
                         |Hdx, Hdy, absVh |                    |
-------------------------|----------------|--------------------|
from (X_accel , Y_accelfrom (X_accel , Y_accel) |Addx, Addy      | X_accel1, Y_accel1 |
                         |Adx, Ady, absVa |                    |
-------------------------|----------------|--------------------|
 

load(fname);



function [prm_mean] = get_prm_day_mean(paramnum,result);
USE: function [prm_mean] = get_prm_day_mean(param_number);
Sharon and Yoram Feb 99

function growsegment(mov, iterations, startp)
USE: function [error_total, params , startp, endp]=growsegment(mov, iterations, startp)

function halfpwrfit(mov, parnumber)
USE: function [params, err]=halfpwrfit(mov, parnumber)
HALFPWRFIT Performs 2/3 power law fit and error calculation
Uses method suggested by Moshe Abeles: half of the points are
used for fitting and the other half for error calculation.
[params, err]=HALFPWRFIT(mov, parnumber)
The formula being fitted is:
        tangvel=K*(radcurv/(1+alpha*radurv))^betha
params(1)=alpha
params(2)=betha
params(3)=K
The second parameter is the number of parameters kept fe number of parameters kept free
during the fit.
parnumber=1: alpha is set to 0, betha to 1/3 and K is fitted.
parnumber=2: alpha is set to 0.05 and betha and K are fitted.
parnumber=3 (default) all 3 parameters are kept free.
 

AUTHOR: Written by Vladimir Litvak 22/10/98



function histmulti5(y,bins)
USE: function [n,x]=histmulti5(y,bins)

DOES: For multidimensional histograms for Matlab 5.

[n,x,nbins]=histmulti5(y,bins)

Every row in y represent a different data point and the columns
are the different coordinates.

For one-dimensional data you should use a column-vectors.
INPUT:
bins: a rowvector => This is the number of bins in that coordinate.
bins: a matrix => The lower limits for the different bins as produced in x.
   Pad at end of the columns with NaN if different number of bins
   in different coordinates (and unequal bins).
The upper limits are the lower limit for the next bound and inf for
the last one.
 

The default is that bins=[10 ... 10];
OUTPUT:
Result:
  n - a multidimensional vector with counts.
  x - Corresponds to the bins-matrix (with NaN-padding if necessary).

The algorithms are constructed to be very fast, and can be superior
to hist even for one-dimensional histograms.

It uses a mex-file to get good performance, and the first tiood performance, and the first time it
automatically compiles the mex-file. If you have trouble compiling
c-files please compile the c-file by hand or erase it.

Rewritten for Matlab 5 by Hans Olsson, Hans.Olsson@dna.lth.se
Original Matlab 4 routine written by Hans Olsson, Hans.Olsson@dna.lth.se.
Inspired by hist2 written by Kirill K. Pankratov, kirill@plume.mit.edu
Which in turn was based on routines written by Hans.Olsson@dna.lth.se



function reply=isscrible(fileName)
DOES: Reads the file header to determine
whether it contains scrible trial.


function j2pos1(jx , jy)
DOES: Converts jdas reading to handle position
CALLING: [xh , yh] = j2pos (jx , jy)
INPUT:  jx, jy  are the readings of JDAS
OUTPUT : xh, yh  are the position of the handle in cm. (0,0) is
                the place of the encoders
accuracy of conversion +/- 0.2 mm
USES: v2pos.

NOTES: Compensation for Gain factor included:
AUTHOR:   SG, 16.9.97
LAST MODIFIED: 16.9.97
SEE ALSO:  j2v



function j2v (jx , jy)

USE: [vx , vy] = j2v (jx , jy)
INPUT:  jx, jy  are the readings of the JDAS A-to-D converter
       &n;              circuit, in Volts.
OUTPUT:  vx, vy  are the voltages (in Volts) at the input
                of the JDAS

constants:
7     --   reading of JDAS for 0 volts (source SG)
0.0024414 -   10V/4096 = gain of ADC
AUTHOR: MA


function jerk_script(filestr,onlymove);
USE: function [jerk_vals] = jerk_script(filestr,onlymove);
[jerk_vals] = jerk_script(filestr,onlymove);
DOES: Calculate and display mean jerk values..

filestr is a string specifying which files to load.
For example show_loc_freq('O*A.*') will load all Ophelai files
in the directory.
AUTHOR: Feb 99 YBS.

A flag to indicate if data should be resampled so that
equal samples are made in space and not in time.



function longSeq(initialList)
DOES: longSeq  find the longest rub of sequential indices
USE: longestSeq = longSeq(initialList)

AUTHOR: Jan 1999 MA


function make2x2
USE: aH = make2x2
DOES: to generate a two by two figure with axes
filling up most of the plot.
aH  - an array of 4 handles for the axes

AUTHOR: mar1999 MA

function aH = make2x2


function make_targ_ WIDTH="100%">function make_targ_prob(xbins,ybins);
USE: function [probmat] = make_targ_prob(xbins,ybins);
DOES: Generate a matrix of probabilities of target locations.
Accepts the x and y locations of the 'location matrix' - and using the
target centers and radii, calculates the probability of getting a target
being in each place.
AUTHOR: YBS Feb 1999.
Change to make correct normalizations for matrices.

function make_targets_loc(size,y_shift,x_shift,targetrad);
USE: function [target] = make_targets_loc(size,y_shift,x_shift,targetrad);
DOES: RETURN CENTERS OF TARGETS.
function plottargtes = plottargets(size,y_shift,x_shift,targetrad);
plots 19 targets that make up the grid.
INPUT: size - horizontal distance between adjacent targets. This is the parameter
  given in the params.dat file ( of the bc) as HOR_DIST.
y_shift - General offset of the woprkspace. Given in cm. Defind by parameter
TAR_OFFSET in the params.dat file ( of bc).
targetrad - radius of target. This should be given in mm. This appears in the
dat files as target_rad ( the name changes depending on the trial).
The function toggles 'HOLD' to ON.
Plots all targets ceneters as well as a circle with the same size of the targets
around them.


compact = mdeblank (string)
DOES:will eliminate any 'space' characters from string'space' characters from string
AUTHOR:MA Feb-1998


function meanjerk(mov)
USE: function [meanjerk]=meanjerk(mov)
DOES: Calculate mean jerk of given motion:
Based on the formula in the viviani and flash paper
"Journal of experimental psychology 1995 vol 21 no1,32-53"
AUTHOR: FEB 99 YBS

Cutoff frequency of motion derivation
cutoff = 4;
Filter the movement
INPUT: mov = filt(mov , cutoff)
mov=eqspacedist(mov)


partRes= mExtract(results, str)
DOES:    to generate an array with all the i-th values
   of prms in dayMov.
Two forms to this routine: Either both argment are strings,
and then:
 INPUT:
   results - is the dir where the daySeg file is.
     str     - is a string that specifies which column to extract
or results is an array and str is a number,
and then:
OUTPUT:    results - is the name of the array of prms (whos name is
     is usually dayMov).
     str   -  is the index to the column to be used
partRes   -  is a column vector with all the values
In form one, only movements which are contained completely within a file
are considered.  In form two all the parameters in DAYMOV are returned.parameters in DAYMOV are returned.

Mar-1999 from rExtract at 3-Mar-1999.  MA



function monk_sum(parentDir)
DOES: To extract statistical summary of all days.
USE: [summary , sumDescript] = monk_sum(parentDir)
it is assumed that under parentDir there are only
directories with a valid date in them, and that eack of them
contains the DAYSEG and DAYRES data bases.
parentDir   - the dir where all the data is available in subdirectories
summary     - a (Nx45) array with summary data
sumDescript - description what is in each column of data

NOTE this code is inefficient since mExtract (rExtarct) will read
   the data bases again and again...

AUTHOR: Mar-1999 MA



script mov2mat (inFile , outFile)
DOES: predefine  inFile, than run
  >> mov2mat

  To convert cobtinuous analog channels in "Opher" format to
  *.mat files suitable for motion analysis.
  inFile  - full name of the input file (e.g. 'o01dec7a.001').
            Only the data between Start and End will be converted
  outFile - the beginnig of the MAT file (.e.g. 'o01dec7a_001').
  The procedure will generate an output file with the following
  variables (where N is no_of_samples):
HandX &f_samples):
HandX      Nx1  X-position in mm
HandXdt    1x1  sampling interval in secs.
HandY      Nx1  Y-position in mm
HandYdt    1x1  sampling interval in secs.
X_accel    Nx1  acceleration in direction perpendicular to Ulna (m/s^2)
X_acceldt  1x1  sampling interval in secs.
Y_accel    Nx1  acceleration in direction parallel to Ulna (m/s^2)
Y_acceldt  1x1  sampling interval in secs.
alpha      Nx1  angle of shoulder in degrees
alphadt    1x1  sampling interval in secs.
beta       Nx1  angle of elbow in degrees
betadt     1x1  sampling interval in secs.
M.A. Feb-1998
additions
timeList   1x1  Time of first sample relative to start of file
endS       1x1  No. of samples in the record
endData    1x1  Time of end of file
hour       1x8  Time of day for this file
Notes:   END DATA is not put right by Yevgeni's translation.  He
   puts the JDAS clock count instead of the sampling count.  This
   is circumvented in this version. MA.

Uted in this version. MA.

UPDATES
Jan-1999  preallocate arrays according to:
           N=1.1*FileSize/8/NoOfChannels.  M.A.
  ???      No OutFile is generated (Vladimir?)
  ???      Variable names are as defined in the analog data
           input file. (Vladimir?)



function mov2mat (inFile , outFile)
parameters for identifying key-words


function moveOnly (mov, delay)
DOES: moveOnly  to animate only the motion with the rewards
INPUT: mov - the movement object with reward times (if relevant)
delay a parameter that depends on the Hardware and slowes
       down the animation.

this is based on PrepareAnimate and animateMove
AUTHOR: Feb-1999 MA

UPDATES: Added the display of the targets with the function plottargets
YBs feb 99.



script myatan
MYATAN2 to find arc-tang in range [0-2pi)
angle = myAtan2(y,x)
  y = array with ordinate values
  x = array with abscissa values

AUTHOR:Jan 1999 M.A.


function md = myMAD(x)
myMAD median of absolute deviation from median

function new_show_loc_freq(f>function new_show_loc_freq(filestr,onlymove);
USE: function [curslocs] = new_show_loc_freq(filestr,onlymove);
DOES: Plot the cursor location probability distribution.
filestr is a string specifying which files to load.
For example show_loc_freq('O*A.*') will load all Ophelai files
in the directory.
AUTHOR: Feb 99 YBS.

A flag to indicate if data should be resampled so that
equal samples are made in space and not in time.

UPDATES: changes to graphics by MA:
square bins, bin boundery not shown, title includes date

more changes by YBS - Plot target location as lines and not as markers.
All matrices are normalized to be of area = 1.
Calculate ( and plot ) the probability of getting a reward at each location.
Calculate ( and plot ) the probability of getting a reward * probability of being
in a location - which is proportional to the number of rewards in every location.
Also - the sum of this matrix is calculated and displayed as a measure of succes.
 



function old2new(inFile)


function parse_events_file(filename, eventNumber,qualNumber)
USE: function [eventTimes, quals]=parse_events_file(filename, eventNumber,qualNumber)
DOES: PARSE_EVENTS_FILE Reads events (digital) file and returns event times.
PARSE_EVENTS_FILE(filename, eventNumber)
INPUT: filename - name of the file to reilename - name of the file to read.
eventNumber - number (as appears in the file) of the event of interest.

AUTHOR: Written by Vladimir Litvak Aug-Sep 1998.

Limitations:
  1.  it is assumed that events are separated by white spaces
      and within events code,qual,dt are separated by commas.
  2.  No hexadecimal conversions
  4.  No treatment of qualifiers
UPDATES
  Jan-1999
     lines with comments are ignored  MA

Update:
FEB - 1999
Add the possibilty to specify qualifiers in addition to the events.
This version should be downward compatible - that is a result will
Also be given for events without qualifiers.
YBS.
If no qualifier is given - then the qualifier of the event is returned.


function PrintSeg(colortype)

DOES: PrintSeg (COLORTYPE)  to print a segmentation picture on Post-Script
printer
INPUT: COLORTYPE - 'B' for black and white
            'C' for color
select the figure to be plotted before calling


script playAfile
DOES: playAfile  script to play scribling
movements in one file.
  This is used in conjunction with animateMove
  see explanation of needed varaibles there
  WARNING
  re
  WARNING
  if rewards are too close, an error message is issued
  2 method of retracing are given:
    comet   for slow PCs
    retrace for fast PCs
  comment out the appropriate groups of lines to use one

AUTHOR: Jan 1999 M.A.
Updates
   May-1999   A few graphis enhabcememnts.  REWARD no. is not outputed
              since that erases the dots around it.


script playStrk
does:playStrk  script to play strokes and show power law
  This is used in conjunction with animateMove
  see explanation of needed varaibles there

AUTHOR: Jan 1999 M.A.
UPDATES: May-1999  A few graphic enhancements. b,K,r2 are writen on power-low


script playTgv
DOES: playTgv  A script to play the Tangential velocity vs time
  This is used in conjunction with animateMove
  see explanation of needed varaibles there

AUTHOR: Jan-1999  MA
UPDATES:   May-1999  A few graphic enhancements


function plot1file (file, axisH)

DOES: plot1file (file) to plot the scribling of one file
INPUT: file  - if a string, file is loaded and then filtered
       - else file ibsp;   - else file is a name of movement and it is plotted
axisH - an optional parameter for the axis on which to plot



function  plot1mov (dirName, fileNo, movementNo, oldH)
USE: aH = plot1mov (dirName, fileNo, movementNo, oldH)
DOES: to plot 1
   movements in one directory
   the list of all these names is returned in cat
INPUT: dirName  - the name of the directory to plot
fileNo   - if supplied this file is plotted.
            NOTE file no. is its number in daySeg data base.

movementNo - no. of movement within this file
oldH     - Axis Handle where to plot
            if not provided - start a new figure
OUTPUT: aH       - handle to the axis on which plotted
the program assumes that daySeg data base is available in the
same directory

AUTHOR: Mar 1999 M.A.


function plot1movWithTV (dirName, fileNo, movementNo, strokes)
USE: mov = plot1movWithTV (dirName, fileNo, movementNo, strokes) to plot 1
 DOES:   movements in one directory with its tangential velocity
   the list of all these names is returned in cat
INPUT: dirName  - the name of the directory to p the name of the directory to plot
fileNo   - if supplied this file is plotted.
            NOTE file no. is its number in daySeg data base.
movementNo - no. of movement within this file
strokes    - an optional list of stroke numbers.
              if provided a second figure with power-law
              fit is drawn.
mov      - the movement
the program assumes that daySeg data base is available in the
same directory

AUTHOR: Mar 1999 M.A.


function plot1prm(summary, sumDescript, fieldName, prmNo, days)
DOES: plot1prm  to plot 1 parameter from summary of all data bases
USE: handle = plot1prm(summary, sumDescript, fieldName, prmNo, days)
INPUT: summary   - is an array of structures, as generated by MONK_SUM
sumDescript - text array with description of variables
fieldName - a string with the field to plot
prmNo     - no. of parameter within this field (see sumDescript)
style     - string saying how to plot:
           - 'dots'      dots with no connecting line
 &nbsnecting line
           - 'line'      dots with connecting line
           - 'bar'       bar chart
           - 'stem'      stem graph
           - 'errorbar'  connected bar with +/- 2 SE
           if omitted STEM is used
days      - list of indices to plot. If ommited, all are plotted
handle    - a handle to the axes
example:
   plot1prm(summary, sumDescript , 'median' , 3, [1, 2, 3]);

AUTHOR: Mar 1999 MA


script plotAfile
DOES: plotAfile  script to plot the scribling of one file
  This is used in conjunction with animateMove
  see explanation of needed varaibles there

AUTHOR: Jan 1999 M.A.
function plotAllM(fname, inDir)


script plotAllM
DOES: plotAllM  -Plot All Motions :  In the file FNAME
  this is a script that should be used like:
    >>fname = 'ana050';
    >>getManip
    >>toPrint = 0;
    >>figTitle = 'o03aug7p;  >>figTitle = 'o03aug7';
    >>plotAllM

  parameters
  fname is the name of the analog data file (after conversion
        from NDA format to Opher format).
  toPrint 0 if only to plot on screen,
          1 if also to send to printer.
          2 if to save figure as a postscript file
  figTitle  text to be added at the top of each page
  see getManip for the definition of varaibles used here.

  you may modify the 'plotStructure' to accomodate
  different no. of motions per page

  this script uses: get_s_e  function.



script plotTgv
plotTgv  To plot motion and tangential velocity
  This is used in conjunction with animateMove
  see explanation of needed varaibles there

AUTHOR: Jan 1999 M.A.



function plotThresh(mov)

function plot_each_mov(dirName, fileNo, toFileNo)
USE: cat = plot_each_mov (dirName, fileNo, toFileNo) to plot the scribling of all
   movements in one directory (or one file within this
   directory).  Only complete movements are
   plotted.  After each movement the user types its "name".    the list of all these names is returned in cat
INPUT: dirName  - the name of the directory to plot
fileNo   - if supplied this file is plotted.
            NOTE file no. is its number in daySeg data base.
          - if omitted, all files are plotted
toFileNo - if supplied files fileNo-toFileNo are plotted
            if omitted all the files to the end are plotted
cat      - a list of movements categorizations 0-9 simple movements
   two digits or more compound movements
the program assumes that daySeg data base is available in the
same directory

AUTHOR: Mar 1999 M.A.


function plottargetbars(targets,probmat,xbin,ybin);
USE: function plottargetbars = plottargetbars(targets,probmat,xbin,ybin);
DOES: Plot lines for targets to fit for the new2_show_freq_loc
AUTHOR: YBS - Feb 1999


function plottargets(size,y_shift,x_shift,targetrad);
USE: plottargtes = plottargets(size,y_shift,x_shift,targetrad);
DOES: plots 19 targets that make up the grid.
INPUT: size - horizontal distance between adjacent targets. This is the parameter
  given in the params.dat fp; given in the params.dat file ( of the bc) as HOR_DIST.
y_shift - General offset of the woprkspace. Given in cm. Defind by parameter
TAR_OFFSET in the params.dat file ( of bc).
targetrad - radius of target. This should be given in mm. This appears in the
dat files as target_rad ( the name changes depending on the trial).
The function toggles 'HOLD' to ON.
Plots all targets ceneters as well as a circle with the same size of the targets
around them.


function powerlawbreak(filename,timewin);
example: powerlawbreak('U19AUG9A.040.mat',100);
DOES: compute the time points of hand trajectory data where powerlaw is violated.
It must be called from the 'motion' subdir/
INPUT: Filename is the '.mat' file as previously compiled from
       the motion data.
       Timewin is the time window in data points within which
       the power law fit is calculated.
OUTPUT: The program extracts and plots
        (1) the tang velocity in time
        (2) the Beta coefficient of the 1/3 power law fit
        (3) the Error b/w the actual trajectory and 1/3 power law fit
        In a separate figure the histograa separate figure the histograms of all of these are generated.
UPDATES: 24/FEB/2000 Zoltan Nadasdy


script prepareanimate
DOES: PREPAREANIMATE a script to prepare data for ANIMATEMOVE
prepare the data to be animated in  m
e.g. by load o03aug7a.050.mat

this script will prepare the following variables:
x           - array with filtered X axis values
y           - array with filtered Y axis values
t           - array with sampling times
tv          - array with tangential velocities
fileName    - name of analog file with data
theta       - threshold for separating motion from rest
radcurv     - radius of curvature (r/(1+a*r)
rewardTimes - times at which reward was turned on
segFile     - data about the big movements in the file
movFile     - data about each of the big movement
resFile     - data on each stroke within each movement

Jan-1999  M.A.



script printc
PRINTC  to print current figure in color


function purgeRes (fileName)
DpurgeRes (fileName)
DOES: to purge empty files from RESDAY


function resAday(inDir, outFile, outDir)
USE: function [dayRes, resDescript]=resAday(inDir, outFile, outDir)
DOES: RESADAY to analyse all files of one day and put the results in RESSTRUCT
calling: [dayRes, resDescript]=resAday(inDir, outFile, outDir)
INPUT: inDir  directory from which to take the files, default
        is current directory
outFile name of file with results. default is AllResult

outDir  directory into which to put the result file, default
        is current directory
dayRes a structure with file-name and the results from
        segments.
resDescript is a string matrix with description of what's in the dayRes table.
    results are saved avery 10 files (to have something in case of a crash).
  If outFile already exists it is oppened and all the files which are already there
  are skipped.

AUTHOR: Jan-1999  M.A.
UPDATES:
Feb-1999
   segDay file is used to select the files tp be


function  pwrlaw(params, xdata)
USE: v=pwrlaw(params, xdata)


script pwrsegments

function resAfile(inFile%">function resAfile(inFile, inDir)
USE: resStruct=resAfile(inFile, inDir)
DOES: RESAFILE to analyse one file and put the results in RESSTRUCT
USE: resStruct = resAfile(inFile, inDir)
INPUT: inFile name of *.mat file with movement (see experiment2mat
        for conversion of analog files into mat files).
inDir  directory from which to take the file, default
        is current directory
OUTPUT: resStruct a structure with file-name and the results from
        segments.

AUTHOR: Jan-1999  M.A.



function radius( C )
USE: function [r , x] = radius( C )
DOES: RADIUS Computes the radius of curvature for three points on a circle.
[r , x] = radius( C )
C - matrix of 3x2 where the X-coordinates are in the first
     column, and the Y coordinates on the second column.
r - radius of curvature.  If not readily computable the value
     of -1 is returned
x - the <x,y> centers of the circles
     If C contains more then 3 points, [r,x] are the averages
see copybook 17 1:18 for algorithm

AUTHOR: Feb-1998  MA
UPDATES:
  jan-1999
    When almost a streight line this algorithm may fail (a
 &norithm may fail (a
       sreight line fits the data MUCH better).  in these
function [r , x] = radius( C )
RADIUS Computes the radius of curvature for three points on a circle.
[r , x] = radius( C )
C - matrix of 3x2 where the X-coordinates are in the first
     column, and the Y coordinates on the second column.
r - radius of curvature.  If not readily computable the value
     of -1 is returned
x - the <x,y> centers of the circles
     If C contains more then 3 points, [r,x] are the averages
see copybook 17 1:18 for algorithm

Feb-1998  MA
UPDATES
  jan-1999
    When almost a streight line this algorithm may fail (a
       sreight line fits the data MUCH better).  in these
       conditions only the midpoint and two extreme ones
       are used. MA



function partRes= rExtract(results, str)
DOES:    to generate an array with all the i-th values
   of results in dayRes.
Two forms to this routine: Either both argment are strings,
and then:
 INPUT:     results - is the dir where dayRes and daySeg files are.
    re.
     str     - is a string that specifies which column to extract
or results is an array and str is a number,
and then:
    results - is the name of the array of results (whos name is
     is usually dayRes).
     str   -  is the index to the column to be used
OUTPUT: partRes   -  is a column vector with all the values
In form one, when a file starts/ends in a middle of a stroke, these strokes
are ignored.  In form two all the parameters in DAYRES are returned.

AUTHOR: Jan-1999 MA

function partRes= rExtract(results, str)


function rotation(mov,i)
USE: ccw = rotation(mov,i)
DOES: ROTATION returns  1 if rotating Counter Clock Wise
                  -1 if CW
                   0 if radius is > 200 mm.
mov is a movement object
i   is the place around which the rotation is requested.

AUTHOR: Jan 1999 MA
function ccw = rotation(mov,i)



to run a few days for segment a day
below are a few examples
inD='f:/result/02Dec97';
outD=inD;
resAday(inD,'dayRes',outD);

outD='f:/y(inD,'dayRes',outD);

outD='f:/result/05aug97';
inD='G:/ophelia/05aug97/analog';
digD='G:/ophelia/05aug97/digital';
run_all_night(inD,outD,digD);

outD='f:/result/10jul97';
inD='G:/10jul97/gnral';
digD='G:/10jul97';
run_all_night(inD,outD,digD);
segAday(outD,'daySeg',outD);
resAday(outD,'dayRes',outD);

outD='f:/result/19jul97';
inD='G:/ophelia/19jul97/analog';  WHY files have the name O19JUN7a.* ???
digD='G:/ophelia/19jul97/digital';
run_all_night(inD,outD,digD);
segAday(outD,'daySeg',outD);
resAday(outD,'dayres',outD);
outD='f:/result/21sep97';
inD='G:/ophelia/21sep97/analog';
digD='G:/ophelia/21sep97/digital'; %WHY file here have the name o21OCT7z...?
run_all_night(inD,outD,digD);
segAday(outD,'daySeg',outD);
resAday(outD,'dayres',outD);
outD='f:/result/24oct97';
inD='G:/ophelia/24oct97/analog';
digD='G:/ophelia/24oct97/digital'; %WHY file here have the name o23OCT7z...?
run_all_night(inD,outD,digD);
segAday(outD,'daySeg',outD);
resAday(outD,'dayres',outD);



function run_all_night(inDir, outDir, digDir)
  run_all_night(inDir, outDir, digDir)
DOES:  to run over one day and convert all files
  with 'scrible' in the event file into *.mat files
  with movement data.
INDIR is the name of the directory where all the analog
   mot all the analog
   motion files exist, default is current directory.
OUTDIR is the directory into which the MAT files will be written
   if omitted the current directory is used.
DIGDIR is the directory where the event files are stored.
   Default is INDIR..\DIGITAL
   if inDir == -1 no digital files were found - donot look
      for events.
assumptions:
   1. analog file names end with A.
   2. event files are in directory ..\digital
   3. event files have the same name like analog files
      but with Z at the end instead of A
   4. scrible appears in "TITLE(3)=..." token within the event-file
the movement will be written to the OUTDIR
and will have the form FILENAME.EXT.MAT.

UPDATES
  Feb-1999
     Files which exist are not converted again
     after conversion, data bases are also formed.
     If digital file does not contain 0,1,0 and 0,2, the translation
     is not carried (no info about rewards).
  Jan-1999
     digDir is sent also to experiment2mat.
     inDir, digDir and outDir added.  M.A.

parameters



function save_data()
save_data
()


function segAday(inDir, outFile, outDir)
USE: [daySeg, dayMov, segDescript, movDescript]=segAday(inDir, outFile, outDir)
DOES: SEGADAY to analyse all files of one day and put the results in
SEGSTRUCT and MOVSTRUCT
[daySeg, dayMov, segDescript, movDescript]= segAday(inDir, outFile, outDir)
INPUT: inDir  directory from which to take the files, default
        is current directory
outFile name of file with results. default is AllResult

outDir  directory into which to put the result file, default
        is current directory
daySeg a structure with file-name and the list of big movements
        within each file.
*Descript are string matrices with description of what's in the day* tables.
    results are saved avery 10 files (to have something in case of a crash).
  If outFile already exists it is oppened and all the files which are already there
  are skipped.

Jan-1999  M.A.
UPDATES
Feb-1999
  only files with less then 800,000 bytes are processed.



script scorePlot
DOES: scorePlot To plot Scores for getting a rewards of all days.
[scores , dates] = scorePlot (parentDir)
it is assumed that under parentDir thessumed that under parentDir there are only
directories with a valid date in them, and that each of them
contains the dateFREQ.mat file.
INPUT: parentDir   - the dir where all the data is available in subdirectories
summary     - a (Nx45) array with summary data
sumDescript - description what is in each column of data

NOTE this code is inefficient since mExtract (rExtarct) will read
   the data bases again and again...

AUTHOR:Mar-1999 MA



function scribbleFiles=show_trials()


function segFile(mov, visual)
SEGFILE to describe the major motions within a file
[bigMov, movPrm] = segFile(mov, visual)
  to show how mov is parsed into compund motions
  MOV is the name of a movement object
  VISUAL if 'vis' pictures of the motion are generated
  BIGMOV Structure with information on the big motions
   .fileName
   .prms(1)  Threshold for parsing
   .prms( 2) min duration below threshold
   .prms( 3) No. of movements
   .prms( 4) 1-st motion started inside the file
   .prms( 5) last motion ended inside file
   .prms( 6) time of last sample (~file's duration)
   .prms( 7) No. of rewards
   .prms( 8) mean time between rewards
&nbn time between rewards
  MOVPRM structure with parameters for each movement in the file
   .fileName
   .prms(i, 1)  No. of motion
   .prms(i, 2)  start X
   .prms(i, 3)  start Y
   .prms(i, 4)  end X
   .prms(i, 5)  end Y
   .prms(i, 6)  center of "mass" X
   .prms(i, 7)  center of "mass" Y
   .prms(i, 8)  total length (mm)
   .prms(i, 9)  angular length relative to ceneter of mass (radians)
   .prms(i,10)  time of start of motion (in the file)
   .prms(i,11)  time of end of motion (in the file)
   .prms(i,12)  duration (seconds)
   .prms(i,13)  No. of rewards
   .prms(i,14)  Time of 1-st reward (relative to start of motion)
   .prms(i,15)  Time of last reward (relative to start of motion)

WARNING
   mov.file should contain a name like o03aug7...

Jan-1999  M.A.
UPDATES
Feb-1999  Problems with m(1).x when m is only 1x1. --fixed
      the ugly way.  MA



function segmanlz(mov)
USE: result=segmanlz(mov)
DOES: SEGMANLZ Analyzes a segment of movement and returns the results in a vector
result=SEGMANLZ(mov) result=SEGMANLZ(mov)
The result vector contains (in this order): movement starting time
(relative to the beginning of the file), movement duration, movement
length (is space), power law betha, power law K, r^2 for the fit between
log(rstar) and log(tangvel).

AUTHOR: Written by Vladimir Litvak 10/11/98
UPDATES
       General Power Law may be skipped by defining NoGenPwrLaw
  Feb-1999
      mean absolute jerk is added
      mean velocity and mean acceleration are added
      treatment for case that they are not available.
  Jan-1999
      info about direction of rotation at center of segment
           is added (1 means counter clock wise),
      median radius of curvature is added.


function  segments(mov, analsTypes, visual)
USE: result = segments(mov, analsTypes, visual)
DOES: SEGMENTS Perform different analysis on movement segments and
stores the result in an ASCII file.
result=SEGMENTS(mov,  analsTypes)
analsTypes - vector of numbers indicating how to break the movement
           into segments.
1 - from minimum to the next minimum of imum to the next minimum of tangential velocity.
2 - from minimum to the next maximum of tangential velocity.
3 - from maximum to the next minimum of tangential velocity.
4 - between curvature breakpoints.
5 - between points of maximal distance from the center of mass.
6 - both min-to-max and max-to-min (2+3 above).
Outputs:
  RESULT array with description of parameters of strokes. It contains
       1 Date coded in binary
       2 extension no. of file
       3 analsType (from input)
       4 # of movement within file
       5 # of stroke within the movement
       6 start time of stroke
       7 duration of stroke
       8 length of stroke (in mm)
       9 alpha  | all three parameters as
      10 beta   | extracted by generalizad power low
      11 K      | V=K*(R/(1+alpha*R)^beta
      12 beta   with alpha fixed (currently at 0.005 1/mm)
      13 K            "
 &nsp;  "
      14 Rsquare correlation coeficient squared
      15 ccw   -1 if motion is cw, 0 if streight line, 1 if ccw
      16 median radius of curvature
      17 X componnet of mean velocity
      18 Y             "
      19 X component of mean acceleration
      20 Y             "
      21 a split storke - 0 as defined by segments
                          1 further splitted due to bad fit (not implemented)
      22 median of absolute jerk
      23 mad of absolute jerk
      24 value of fixed alpha (for 12,13) fixed at 0.005
WARNING
   mov.file should contain a name like o03aug7...

UPDATES
Apr-1999  if VISUAL 3 user responses are allowed:
    <Enter> continue to next segment
    N       jump to next major movement
     &nnbsp;         Q       quit the program
Feb-1999
     problems when only one segment:
     size(segments.x) = nnn  BUT size(segments(1).x) = 1 --
     by-passed the hard way.
     special test for files with absolutely no motion
Jan-1999
      Using the same figure if exists.
      printout of K supressed
      analsType 6 was added  MA
      If most (80%) of the file is without motion - ignore the
      file (is Sleeping)   MA


To break daySeg and dayMov to simple table
  to build a very long tables with 1:1 match among indices
[fmtable, mtable, fstable, stable] = seg_table(daySeg, dayMov)
daySeg, dayMov - optional names where the segments and movements
    of one day are stored, if either is not supplied the DAYSEG
    data-base is loaded from the current directory
fstable - a list of the file names in daySeg
stable  - a table with the daySeg.prms data
mtable  - a table with the dayMov.prms data
fmtable - a table with file-names expanded to have as
    many rows as in mbsp;  many rows as in mtable.

Feb-1999 MA
   WARNING data which is not available is entered as NaN



function seg_table_c_only(daySeg, dayMov)
DOES:  To break daySeg and dayMov to simple tables
  to build a very long tables with 1:1 match among indices
  only those files (and data) which contains full movements
  are included.
USE: [fmtable, mtable, fstable, stable] = seg_table_c_only(daySeg, dayMov)
INPUT: daySeg, dayMov - optional names where the segments and movements
    of one day are stored, if either is not supplied the DAYSEG
    data-base is loaded from the current directory
fstable - a list of the file names in daySeg with at least one complete
           movement.
stable  - a table with the daySeg.prms data. prms(3,4,5) are adjusted
           for complete movements only.
mtable  - a table with the dayMov.prms data. Only data about complete
           movements are included.
OUTPUT: fmtable - a table with file-names expanded to have as
           many rows as in mtable.
   NOTE in stable only params 3,4,5 are modified
 &,5 are modified
   WARNING data which is not available is entered as NaN

AUTHOR: Feb-1999 MA



function show_loc_freq(filestr,onlymove);
USE: [curslocs] = show_loc_freq(filestr,onlymove);
DOES:Plot the cursor location probability distribution.
filestr is a string specifying which files to load.
For example show_loc_freq('O*A.*') will load all Ophelai files
in the directory.
AUTHOR: Feb 99 YBS.

A flag to indicate if data should be resampled so that
equal samples are made in space and not in time.
resample= 1;


function show_xy_freq2(filestr,onlymove,resample);
USE: [curslocs,real_target_mat,reward_prob,score,x,y,run_error] = DOES: show_xy_freq2(filestr,onlymove,resample);
[curslocs,real_target_mat,reward_prob,score] = show_xy_freq2(filestr,onlymove,resample);
Plot the cursor location probability distribution.
 

INPUTS:
filestr is a string specifying which files to load.
For example show_loc_freq('O*A.*') will load all Ophelai files
in the directory.
onlymove - Take only movements (default 1).
resample - A flag to indicate if data should be resampled so that (default 0)
equal samples are made in space and not in time.

OUTPUTS:
curslocs         - matrix with probabilities of location of cursors
real_target_mat  sors
real_target_mat  - the matrix of the probabilities where the targets appears.
reward_prob      - the bin by bin product of both matrices.
score            - the sum of the reward_prob matrix
x,y - the x and y bins of the matrices.
Error is 1 if program aborted because of no data.

Mar 1999
Error: when hitting a file with no motion the process is stopped
     FIXED by skipping files where no motion occures M.A.
     all outputs are initialized.
Feb 99 YBS.

changes to graphics by MA:
square bins, bin boundery not shown, title includes date

more changes by YBS - Plot target location as lines and not as markers.
All matrices are normalized to be of area = 1.
Calculate ( and plot ) the probability of getting a reward at each location.
Calculate ( and plot ) the probability of getting a reward * probability of being
in a location - which is proportional to the number of rewards in every location.
Also - the sum of this matrix is calculated and displayed as a measure of succes.

More changes: 21/2/99 YBS
Currect normalization of matrices:



function sinwave(params, t)
USE:x=sinwave(params, t)


function slidewindow(mov, winsize)
USE:  [betha, K, err]=s)
USE:  [betha, K, err]=slidewindow(mov, winsize)function segments=split(mov, startPoints, endPoints)
DOES: SPLIT Splits the movement into segments.
segments=SPLIT(mov, startPoints, endPoints)
Returns an array of movement objects each of which is a segment
of the original movement.
INPUT: mov - movement to split
startPoints - starting points of the segments. If endPoints argument is omitted
               the movement is divided into segments between these points (including
               the beginning and the end.
endPoints-    ending points of the segments. The numbers of starting and ending points
               should be equal.

AUTHOR:  Written by Vladimir Litvak 1/10/98
   WARNING  only fields with column-length
          equal to that of mov.npoints are being split

UPDATES
  Feb-1999
     problems when only one segment:
     size(segments.x) = nnn  BUT size(segments(1).x) = 1 --
     by-passed the hard way.
     since subrefs does not treat 2-dim arrays, these


function split(mov, startPoints, endPoints)
USE: segments=split(mov, startPoints, endPoints)
DOES: SPLIT Splits the movement into segments.
segments=SPLIT(mov, startPoints, endPoints)
Returns an array of movement objects each of which is a segment
of the original movement.
INPUT: mov - movement to split
startPoints - starting points of the segments. If endPoints argument is omitted
               the movement is divided into segments between these points (including
               the beginning and the end.
endPoints-    ending points of the segments. The numbers of starting and ending points
               should be equal.

AUTHOR:  Written by Vladimir Litvak 1/10/98
   WARNING  only fields with column-length
          equal to that of mov.npoints are being split

UPDATES
  Feb-1999
     problems when only one segment:
     size(segments.x) = nnn  BUT size(segments(1).x) = 1 --
     by-passed the hard way.
     since subrefs does ;  since subrefs does not treat 2-dim arrays, these cases are treated
     specificaly here.  MA
     FILE and firstSample are treated here individualy
     EVENTS are copied with the ones within the range only
     if EVENTS has more then two cells the other are copied too
     only for the time range which falls in the segment.
cases are treated
     specificaly here.  MA
     FILE and firstSample are treated here individualy
     EVENTS are copied with the ones within the range only
     if EVENTS has more then two cells the other are copied too
     only for the time range which falls in the segment.



function v2pos(vx , vy)
DOES: convert voltages to position of handle
USE: [xh , yh] = v2pos (vx , vy)
INPUT:  vx, vy  are the voltages as read from the encoders
                 circuit, in Volts.
OUTPUT:  xh, yh  are the position of the handle in cm. (0,0) is
                 the place of the encoders
accuracy of conversion +/- 0br>accuracy of conversion +/- 0.2 m
constants:
  0.314159 - 2*pi*4096/8000/10.24
  26.6     - length of upper arm
  35.9     - length of lower arm
1.5677 , 3.377 are neede to bring the position to the calibration point
             when vx=1.99 , vy=1.275
function viewer3d(command,parameter1)
  Template for uicontrol form.

  If you rename this file:
    1. Replace all occurrences of 'viewer3d' with desired 'filename'.
    2. Resave this file as that 'filename'.

  To edit callbacks :
          Fill in the 'elseif command == ...' sections with
          code that does something other than what was
          defined in GUIMaker 2.0. Be creative!
          The code must lie between the CALLBACK_START and
          CALLBACK_STOP comments for each objects callback
          command section if you want to be able to edit
          the GUI with GUIMaker later on and not lose the
  &t lose the
          callback information.


 function viewer3d(command,parameter1)
  viewer3d.m

USE:  viewer3d(command,parameter1)

           where

               'command' is used to perform some defined operation
               'parameter1' can be used pass extra information
                              to a defined operation

           If no arguements are used the GUI is initialized.
 

  Created: 8-May-95
  Using  : GUI Maker Ver 2.0 by Patrick Marchand
                         (pmarchan@motown.ge.com)
  AUTHOR:
  UPDATES:

  Copyright (c)
       Permission is granted to modify and re-distribute this
code in any manner as long as this notice is preserved.
All standard disclaimers apply.


<


 
 
/html>