        EMF11P DEVELOPMENT NOTES

        File:   EMF11P.TXT
        Author: Gregory A. Miller
        Date:   4/2/88

------------------------------------------------------------------------

            OVERVIEW/SUMMARY

        This file contains notes on the conversion of a version of Gabriele
Gratton's eye movement correction procedure (EMCP), developed in Emanuel
Donchin's Cognitive Psychophysiology Laboratory (CPL), to run under DEC
Fortran IV. (While I run it under RT-11 and TSX, it should run without
modification under RSX or RSTS/E Fortran IV.) This text file is not well
structured as a programmer's guide (see EMF11P.GID for that). It is better
seen as a historical record of the conversion for users who wish to understand
(and second-guess) its particulars. Nothing here is conceptually difficult,
but there is much material. Included are some comments on how EMCP works, as I
have come to understand it from experience with the algorithm and, in some
cases, from discussions with Gratton. Consultation with him during the
conversion has been invaluable, but I bear full responsibility for decisions
and possible errors or inconveniences in the conversion.
        I first adapted Gratton's program for my own use on the Harris H800
computer in the CPL (F66 compiler, essentially a Fortran IV enhanced as a
structured language). I worked from the 5/13/83 version of Gratton's EMFAST
program, one of several versions he had developed for the Harris. There were
minor changes and enhancements for four different data sets for which I
adapted his program, with my Harris programs named EMF1, EMF2, EMF3, and EMF4
(differing, e.g., in array dimensions and some details of parameters). Each
version was validated against the previous version as changes were introduced.
        The result of the conversion of EMF4 to DEC Fortran was called
EMF11A.FOR. That source file documents in detail the changes made to the
Harris code, and a copy of that source file can be made available if you care
about those details. EMF11A was validated directly against EMF4 on the Harris.
However, EMF11A is not pretty -- it is a hash of two different programmer's
coding and commenting styles, and it preserves a number of features likely to
be used only within the CPL. The only reason a subsequent user would want to
develop his or her own EMCP program on the basis of EMF11A is that, unlike its
progeny, it does not use the '!' feature of DEC Fortran for inserting comments
on instruction lines. Thus, there would be some small convenience in
converting the EMF11A source code to a non-DEC compiler.
        EMF11A begat EMF11B to deal with another data set under RT-11/TSX.
Whereas a main principle in the creation of EMF11A from EMFAST was to minimize
coding changes, I made a number of non-significant changes from EMF11A to
EMF11B. All but the most trivial are documented here, below. In several ways,
EMF11B remained less general than it might be. It was not convenient to
validate it directly against the Harris EMFAST program. EMF11B was validated
against EMF11A.
        EMF11B begat EMF11C. EMF11C is less idiosyncratic and includes the new
option of facilitating correction for HEM instead of VEM artifact.
Essentially, the user may turn off blink detection and correction, leaving the
saccade correction code intact. In my lab, when I collect both VEM and HEM and
have non-midline EEG channels, I pass each subject's data through EMF11C
twice, first to remove VEM from all channels (including the HEM channel), then
to remove HEM from the EEG channels. This method is somewhat similar to what
Gratton developed independently and is based on our discussions of the
problem. We have not compared the two programs directly. EMF11C in
VEM-correction mode was validated against EMF11B. EMF11C in HEM-correction
mode was validated by study of pre- and post-correction plots. Anyone adapting
the program for their own use will want to spend time making such empirical
comparisons.
        EMF11C begat EMF11D, motivated by yet another data set and made more
general. EMF11P is a version of EMF11D polished up for public consumption.
Because EMF11C is the least idiosyncratic and the most flexible, I recommend
it rather than Gratton's EMFAST or my earlier programs as a basis for someone
else's adaptation of the algorithm for their own uses.
        It's impossible to say what changes another user will want to make.
The first step, if a DEC Fortran compiler is available, should probably be to
get the source code to compile and link as is, before any changes are made.
Beyond that, it seems inevitable that another user will want to change array
sizes (e.g., # points per channel per trial), balance of virtual and
non-virtual arrays (if so, be sure to adjust CLEAR and VCLEAR calls
accordingly), I/O format, and file-naming conventions. Users are likely to
need other customizations and to discover optimizations neither Gratton nor I
perceived or implemented.
        The remainder of the present file is a series of documentation notes
on the lineage of the software. You should choose the sections worth studying,
depending on your purpose.
------------------------------------------------------------------------------

>>>>>>  EMF11A notes:

  5/13/85  G. Gratton
        Date of EMFAST program version from which initial RT-11 adaptation
  was done. EMFAST implements a version of the EMCP method described in
  Gratton, Coles, & Donchin (1983, EEG journal).  Selected comments from
  that Gratton program are reproduced here:

              UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN
                      DEPARTMENT OF PSYCHOLOGY
                 COGNITIVE PSYCHOPHYSIOLOGY LABORATORY
                           EMFAST program
                       AUTHOR: GABRIELE GRATTON
                         DATE: MAY 13TH, 1983
             This program is a version of the EMCP program
             that corrects for vertical eye-movements artifacts.
             A description of the procedure can be found in
             Gratton, Coles, and Donchin (1983).  The
             program is a special version that uses
             computational formulae.

  9/16/85  G. Miller
        Created EMF11A -- adaptation of Gabriele Gratton's Harris F66
  Fortran program EMFAST (5/13/83 version) to DEC Fortran IV.
  Note that that Gratton version did not include the ability to recover
  (estimate) off-scale data under some circumstances -- Gratton added
  that capability to his software later than this version.  This DEC
  version doesn't do that recovery, either.

Compile via:  FORT/EXTEND/UNITS:8 EMF11A
Link via:     LINK EMF11A  or  LINK EMF11A,VIRP

------------------------------------------------------------------------------

>>>>>>  Main changes in conversion of Harris EMFAST to DEC EMF11A:

-- Harris structured Fortran (FOR/END FOR & IF/THEN/END IF)
   replaced by traditional Fortran DO's & IF's.  This involved
   adding many line labels.  All labels added are between 5000
   and 5999.  (Though note that line numbers in subroutines which
   are entirely new to the DEC version are not limited to this
   range of line labels.)
-- Array declarations have been modified to declare some arrays
   VIRTUAL to cope with memory limitations.  The Harris Fortran
   SPECIAL COMMON declaration (which accomplishes something like
   virtual arrays in DEC Fortran) has been deleted.
-- CALL VCLEAR has been substituted for CALL CLEAR when the array
   to be zeroed is virtual.  CLEAR and VCLEAR subroutines are now
   contained within the present source file, whereas previously
   they were obtained from the Harris subroutine library.
-- Harris library subroutines SELECT and CHECK are replaced by
   local subroutines with the same name and function.  The new
   versions are simpler and less flexible than the Harris code.
   See comments within the subroutines for documentation.
-- Subroutine OPNFIL has been added to open files.  File LUNs
   have been made variables, set by DATA statement in EMF11A.
-- In the conversion, as little Gratton/Harris code was changed
   as possible.  Potential execution efficiencies and comment
   enhancements have been ignored.  However, some trivial changes
   in FORMAT statements are noted explicitly.  As in Harris version,
   no 'IMPLICIT INTEGER (A-Z)' is used.

------------------------------------------------------------------------------

>>>>>>  Modifications to EMF11A subsequent to the initial conversion:

-- Added variable INPOPT, controlling whether input file is '.RR1'
   and unformatted (faster); or '.T1L' and formatted (more general).
   (INOPT was not retained in EMF11B.)
-- Renamed NOPT to OUTOPT (INTEGER*2), for better mnemonic parallel
   with INPOPT.
-- Renamed parameter NDR to MSEC, better reflecting its function.
-- Made consistent the use of ID(NIDS-4) for NEOG parameter value and
   ID(NIDS-3) for ADCRIT parameter value.

------------------------------------------------------------------------------

>>>>>>  Comments on EMF11A:

        The Harris version of EMCP has been in use for some years in the CPL,
so much experience in its use has accrued. This DEC version is relatively new;
despite its producing the same output as the Harris version, in the cases
evaluated, it is possible that flaws or idiosyncracies exist in it. It should
be used with caution.
        The #1 principle in the conversion from Harris to DEC Fortran was
minimization of changes to Gratton's code and commentary. Future users will
undoubtedly prefer additional customizations or rationalizations for
convenience. (A few idiosyncracies occur because of this. For example, 'TAB'
characters appear at some points in the file, which is convenient using the
KED editor, but, since the Harris editor converts TABs to spaces, most
oppportunities for TABs are actually just spaces.)
        The minimum needed for conversion consisted primarily of changing the
structured Fortran features available in Harris F66 Fortran into conventional
Fortran IV DO-loops and IF statements. All line labels in the range 5000-5999
were introduced in the conversion. Additional major changes were mainly doing
LUN file assignments within the program, rather than prior to it, and
replacing the CPL/Harris standard SELECT/CHECK trial-classification method
with a simpler, though less flexible, method. Several additional changes to
Gratton's Harris EMFAST code had been introduced in Harris versions of the
program (EMF1, EMF2, EMF3, & EMF4) prior to the DEC conversion, for
convenience in handling data from my lab. A few more such changes were
introduced during the conversion to RT-11. These should be of fairly general
use. All lower-case characters were introduced during the conversion; however,
some changes were made in upper case to avoid odd-looking code.
        One change worth noting here is from parameter LAB to parameter
ADCRIT. In Gratton's EMFAST code, LAB was a switch with which the user
specified which of two CPL labs the data were collected in and hence which A/D
converter was used. This choice determined whether 1023 or 2047 would be used
as the criterion for determining whether a data point was off-channel. This
parameter usage made sense given that the EMCP program was used only with raw
A/D data. A more general parameter was needed, however, if one wanted to be
able to process data already converted to uvolts. EMF11A renamed this
parameter to ADCRIT and generalized it, interpreting it as either (a) the
direct edge-of-channel (off-scale) data value to be used for all channels
(e.g., 2047) or (b) a switch indicating that the edge-of-channel value for
each channel was to be obtained from an ID value stored with that channel. The
latter allows for different channels to have different values and for those
values to be different for different trials.
        After initiating the program, the user is prompted to enter two
things: a 6-character subject identifier, which will form part of the file
names of nearly all files used (the full 6 characters must be supplied); and
an indication of whether the input file extension is '.RR1', signifying that
it is an unformatted, sequential-access file or '.T1L', signifying a
formatted, sequential-access file. (The user will probably want to edit these
according to local conventions.) (File-naming conventions in EMF11B and EMF11C
are simpler than this. See below.)

------------------------------------------------------------------------------

>>>>>> EMF11A's assumptions about the structure of the input data file:

(1) Each record holds the data for one channel for one trial.
(2) Each record consists of identifying information (IDs) followed by data
    points.
(3) Records for all the channels of one trial are contiguous.
(4) VEM is the first channel.
(5) EEG channels follow immediately after VEM.
(6) Additional (non-EEG) channels are permitted to follow the EEG channels
    (such channels are passed on to the output file(s) unmodified).
(7) If parameter NEOG is zero, its value is obtained from ID(NIDS-4) within
    the VEM channel record. (This gain ratio between EEG and VEM is used for
    deciding how big an apparent blink needs to be. Note that a single value
    is used for all EEG channels. Thus, the channels need to have been
    recorded at the same gains.)
(8) If parameter ADCRIT is zero, its value is obtained from ID(NIDS-3) within
    each record. (If the absolute value of data point is greater than this
    value, the trial is assumed to be off-scale.)

------------------------------------------------------------------------------

>>>>>> EMF11A output file structure includes:

(1) No VEM channel (the VEM channel, after 'correction', is flat; however, it
    can be included in the output by changing appropriate DO-loops to begin
    with channel 1 instead of channel 2).
(2) Three parameters created by the program are inserted into the last three
    IDs in the averages EMCP writes (unless the IDALL parameter indicates that
    no IDs are to be retained, in which case the three parameters are simply
    prefixed to the data points). These three are bin #, channel #, & # trials
    in average.

------------------------------------------------------------------------------

>>>>>> EMF11A operational notes:

(1) If one data point in one channel is off-scale, or if one 'flat spot' is
    found in one channel, the trial is rejected in entirety.
(2) A 'flat spot' is defined as five consecutive points having the same value.
    This can occur without the data being detected as 'off-scale', for
    example, if the entire trial was actually off scale but prior to EMCP has
    been deviated from baseline, such that all data values are now zero.
(3) The .E3 corrected single-trial output file is opened before other output
    files, to minimize the chance that fragmentation of disk space will cause
    problems. However, a user who could predict the maximum number of blocks
    needed by each output file could declare this explicitly in the OPEN
    statements in subroutine 'OPNFIL'.
(4) Program speed depends on number-crunching power, proportion of array space
    which is virtual, and disk I/O. The latter could be sped up by using an
    unformatted input file (currently an option) and unformatted output files
    (at least the .E1, .E2, and .E3 files; not currently implemented). The
    heavy user should consider converting to entirely unformatted I/O for data
    files. (Not only is the I/O inherently faster, but the files are smaller.)
(5) To facilitate generality, all LUNs have been made variables (file LUN
    values are those assumed by the Harris EMFAST program, however). The user
    should consider whether his/her system is comfortable with the values
    used.

------------------------------------------------------------------------------

>>>>>> EMF11A file usage:

        FilenameLUN (logical unit #)

   Input:       OUT:zzzzzz.T1L  lunin   data
        DK0:EMF11A.SPC  lunsp   specifications

   Output:      OUT:zzzzzz.E1   lune1   corrected averages
        OUT:zzzzzz.E2   lune2   uncorrected averages
        OUT:zzzzzz.E3   lune4   log of run
        OUT:zzzzzz.E4   lune3   corrected single trials

   User:lunti   user terminal input
        lunto   user terminal output

 Waveform output is point w/in channel w/in bin
------------------------------------------------------------------------------

>>>>>>  EMF11B notes:

Compile via:  FORT/EXTEND/UNITS:7 EMF11B
Link via:     LINK EMF11B  or  LINK EMF11B,VIRP

11/7/86  G. Miller
        Created EMF11B from EMF11A -- customization for use with TR1 data. The
.LST file shows 23,040. words of virtual array usage. The .MAP file shows
25,576. words in the .SAV file. (This does not include buffer space needed for
file I/O. Note that a program this size would not run under most RT-11 systems
as is -- more arrays would have to made virtual for RT-11 execution in most
cases.) Under TSX-Plus, the running task is 108KB. All virtual arrays in
subroutines are declared in the main. No local subroutines call other local
subroutines. Note that all lower case, as well as some upper case, added by
Miller. Note that polarity of channels doesn't matter; needn't all be the
same. Much less attempt is made in EMF11B than in EMF11A to be consistent with
Gratton's coding style, though much of his remains.

------------------------------------------------------------------------------

>>>>>>  Main changes from EMF11A to EMF11B:

1)  Array dimensions changed to handle:
        4000    trials
        7       channels (VEM + 6 EEG)
        30      IDs
        100     points/waveform
        5       bins (4 conditions + garbage bin)
2)  Input data file is now INP:zzzzzz.RE instead of OUT:zzzzzz.RR1 or
    OUT:zzzzzz.T1L.  Input format option has been removed.  Input file
    is assumed to be unformatted direct-access (in subroutines OPNFIL and
    STIN, the latter formerly called XREAD).
3)  Input data file is assumed to be INTEGER*2 A/D values, not REAL*4 uvolt
    values.  Corrected single trials are written in same format as single-
    trial input.  Uncorrected and corrected averages are output as REAL*4
    A/D values.
4)  Output files match input files in terms of file type and record size.
    (Note that direct-access file OPEN statements are very particular about
    record size.  Use of options such as "don't retain input IDs" may
    necessitate a different record size.  Run-time computation of correct
    output record size could be automated but hasn't been implemented yet.)
5)  Subroutine STIN (formerly called XREAD) modified (line labels 6000+) to:
        (a) optionally swap channels 1 & 2 (HEM, VEM --> VEM, HEM),
            depending on new parameter read from spec file, 'HEM'
        (b) read INTEGER*2, not REAL*4, but copy into REAL*4
6)  Fewer arrays virtual (& changed relevant VCLEAR calls to CLEAR).
7)  Array XX now INTEGER*2 instead of REAL*4, to save space.
8)  Numerous cosmetic changes in commenting, spacing, etc.
9)  Removal of most indications of old differences between Harris and DEC
    code, to improve readability.
10) NOUT array now data type BYTE, to save space.
11) Array 'X' in STIN (formerly called XREAD) was really 'V' in the main and
    elsewhere in Gratton's code.  This was preserved in EMF11A but is changed
    in EMF11B -- now in XREAD it's called 'V' as elsewhere.
12) Array 'X2' in subroutine AVOUT (formerly called XWR) was really array
    'AVG' in the main and elsewhere in Gratton's code.  This was preserved in
    EMF11A but is changed in EMF11B -- now in AVOUT it's called 'AVG' as
    elsewhere.
13) Changed some array names within subroutines XSUM & CORFAC (formerly called
    XWRITE) to clarify that their locally named arrays refer to different
    arrays in the caller at different times.
14) Changed array name 'BS' in BASELN (formerly called XBAS) to 'BSLN' to
    match caller usage.
15) Removed array 'INDEX' -- purpose accomplished with implied DO.
16) Specification file no longer kept open throughout run, to reduce by one
    the number of files simultaneously open.
17) Subroutine names changed to be more mnenomic.
18) Program queries user for which trials to drop.  NOUT array stores this
    information.  Note that program accesses trials by sequence, not on the
    basis of some internal trial #.  Thus, if the Nth trial is dropped in the
    course of the run, the (N+1)th trial in the input file is positioned as
    the Nth trial in the single-trial output file.
19) Program inputs the 'maxtrl' trials from a series of input files.  Most of
    the code for this is in subroutine STIN, driven by the current record
    number ('nextin') because EMF11B assumes that there is a different input
    file for each trial block.  If, instead, there were just one input file,
    setting up and opening this file could be done within subroutine OPNFIL.
20) Changed most DO loop indices from I-J-K to mnemonic indices (np for NPTS,
    ni for NIDS, nc for NCHAN, nb for NBIN).
21) "Correction factor" for channel 1 (VEM) no longer printed.  The value
    formerly printed there was not a correction factor but was the EEG/VEM
    ratio.
22) System time and date are read and output at various points, to
    facilitate monitoring execution time.
23) Output to log file is abbreviated:  (a) Calls to subroutine TABLE are
    commented out; (b) Listing of rejected trials formatted very differently
    (done by new subroutine OUTLST).
24) Subroutine CORFAC (formerly called XWRITE) now computes & prints actual
    correction factors (regression betas), not factors divided by EEG/VEM gain
    ratios (no change in adjustment algorithm, just a change in what user is
    told in log file).  This permits deletion of the BXF and SXF arrays.
25) Parameter NEOG, which was the EEG/VEM gain ratio, has been
    replaced by EOGSEN, the digitization sensitivity for the VEM
    channel in A/D units per uvolt.  Gratton's EMFAST program on the
    Harris computer and Miller's EMF11A program under DEC RT-11 used
    NEOG, a gain ratio, as a means to determine the minimum size of
    a blink in the VEM channel.  This usage assumed a fairly stable
    (over time) and consistent (over channels) gain in EEG channels.
    Thus, if VEM gain varied, this would be well reflected in a ratio
    of EEG gain to VEM gain.
        However, EEG channel gain per se is of no importance to the
    EMCP algorithm, which simply computes the correlation (and
    resulting regression beta weight) between the VEM channel and each
    EEG channel.  A more direct parameter usage for NEOG would be to
    have it represent VEM channel gain; specifically, how large a
    value in the VEM channel data scanned by the EMCP program should
    constitute a blink.
        The blink-detection algorithm (see subroutine BLINK, formerly
    called XFIND), checks for a slope, not an amplitude, in judging
    the presence of a blink.  So, the NEOG parameter can't simply
    carry the amplitude value of a blink.  The parameter has been
    renamed here to EOGSEN, conveying that it is the VEM-channel
    digitization sensitivity, and is used in a slightly different way.
    This change allows the parameter to be more general.  The issue is
    explained in the following paragraphs.
        If the data presented to EMF11B are in the original A/D units,
    then the units of EOGSEN should be A/D units per uvolt.  If,
    instead, the data presented to EMF11B are in some other units,
    then EOGSEN should be converted in parallel prior to use in
    EMF11B.  For example, if VEM data in the input file have already
    been converted from A/D units to uvolts (by multiplying each data
    point by some conversion factor denominated in "uvolts per A/D
    unit" units), then EOGSEN (in A/D units per uvolt) should also
    have been multiplied by this same conversion factor.  (The result
    is simply a scalar value, without units.)
        Exactly what function of EOGSEN the BLINK subroutine should
    use to judge blink slope is an empirical question which each
    investigator may want to evaluate for each laboratory and, perhaps,
    for each paradigm and for each subject.  In the CPL, where Gratton
    works, EEG channel gain is kept at 10,000, and data are submitted
    to EMCP in raw A/D form.  With NEOG being the EEG/VEM gain ratio,
    he found that the function ICRIT=150/NEOG worked well across
    paradigms and subjects (ICRIT is the slope criterion for judging
    the presence of a blink).  Given a constant EEG gain, changes in
    NEOG reflect changes in VEM gain; decreases in VEM gain produce
    increases in the value of the NEOG ratio and decreases in the
    value of ICRIT.  Thus, ICRIT is directly proportional to VEM
    gain, as it should be -- higher VEM gain would mean, in A/D units,
    a steeper rise/fall (slope) in the data during a blink.
        Instead of ICRIT=150/NEOG, the present program uses
    ICRIT=14*EOGSEN.  The value '14' was chosen simply because it's
    arithmetically equivalent to Gratton's formula for one data set
    on which this was developed.
        The physiological blink, of course, occurs in uvolts, not A/D
    units.  Thus, if the VEM data presented to EMCP have already been
    converted to uvolts, VEM amplifier gain has already been accounted
    for, and NEOG (or EOGSEN) should not vary as a function of gain.
    That's why the EOGSEN parameter must be treated to the same
    arithmetic as the VEM data, which essentially removes the units
    from EOGSEN.
        One final wrinkle:  since EMCP's blink-detection algorithm
    actually uses slope (amplitude change per msec), not simple
    amplitude, the units of the ICRIT slope criterion must be in
    amplitude change per sec.  Since, in Gratton's function
    ICRIT=150/NEOG, the NEOG has no units (gain divided by gain
    cancels to nothing), the '150' has implicit units of some sort of
    amplitude change per msec, with the units of that amplitude, for
    Gratton, in A/D units.  If the VEM data, however, are no longer in
    A/D units -- in uvolts, perhaps --, then the EOGSEN parameter must
    carry that change in units.  (Note:  the parameter MSEC is needed
    for just this reason -- determination of change per unit time, when
    what is available is only change per data point.)
        Given these considerations, in order to generalize the use of
    NEOG it has here been re-conceptualized, not as gain ratio, but
    as VEM channel digitization sensitivity.  (Again, EEG channel gain
    is irrelevant to EMCP except when, in Gratton's EMFAST program,
    it is assumed to be a stable value on which to peg VEM gain.)  If
    the VEM channel is in raw A/D units, then EOGSEN should be the A/D
    units per uvolt sensitivity at which the VEM was digitized; if the
    VEM data have been subjected to some multiplier, EOGSEN must be
    similarly treated.
        Three options exist for specifying EOGSEN:
    1) If the value of EOGSEN in the specification file is greater
    than zero, then that value is used for all trials.
    2) If the value of EOGSEN in the specification file is equal to
    zero, then the REAL*4 value stored in INTEGER*2 IDs 29 & 30 in the
    VEM channel is used as the EOGSEN value.
    3) If the value of EOGSEN in the specification file is less than
    zero, then the inverse of the REAL*4 value stored in INTEGER*2 IDs
    29 & 30 in the VEM channel is used as the EOGSEN value (e.g., if
    the inverse of sensitivity was stored in the IDs -- uvolts per A/D
    unit).
        Options 2 and 3 allow different trials to have different VEM
    digitization sensitivities.
        These options reflect local conventions and convenience; another
    users may want to alter them.
        Note that the value of the EOGSEN parameter is irrelevant when
    the HEM parameter is set to -1 (i.e., no blink detection/correction).
26) Trials on which a blink was detected are now listed (flag stored
    in array NOUT).  The full set of values for what is stored in NOUT:
    -1=OK w/ blink, 0=OK w/o blink, 1=off-scale, 2=flat EEG, 3=dropped by user
27) Replaced subroutine WST -- moved code to main.
28) Deleted unnecessary CONTINUE statements and line numbers.
29) New parameter in spec file, 'HEM'.  The present program corrects
    all channels for VEM artifact, doing nothing about HEM artifact.
    If an HEM channel is included in the data set, HEM data are
    corrected for VEM artifact (correlation) exactly as if HEM were
    an EEG channel.  The 'HEM' parameter controls this:
        0 = no HEM in data set
        1 = HEM is the second channel; exclude the HEM channel from
            flat-spot checking in subroutine BADTRL
        2 = HEM is the first channel; swap it with the second channel
            in subroutine STIN so that VEM is first and HEM is second;
            exclude the HEM channel from flat-spot checking in
            subroutine BADTRL
    Note that, as in Gratton's program, the data set can contain still
    other channels, if these are positioned after the last EEG channel
    -- this is controlled by parameter 'NCHAN1'.
30) Program-generated parameters were inserted into the last three IDs for
    computed averages -- changed from bin #, channel #, & # trials in
    average.  If the user's IDALL parameter indicates that no input IDs
    are to be output, these three program-generated IDs are output, prefixed
    to the data points in the average (unchanged from EMFAST and EMF11A).
    However, in EMF11B, if IDALL indicates that all input IDs are to retained,
    then the only program-generated parameters which are output are bin # and
    # trials in average.  These are inserted in ID(NIDS-3) and ID(NIDS-2).

------------------------------------------------------------------------------

>>>>>>  EMF11B file usage:

     Input:     INP:zzzzzb.RE   lunin   data for block # 'b'
        DK:EMF11B.SPC   lunsp   specifications (ASCII)
     Output:    OUT:zzzzzA.E1   lune1   corrected averages across All blocks
        OUT:zzzzzA.E2   lune2   uncorrected averages
        OUT:zzzzzA.E3   lune3   corrected single trials
        OUT:zzzzzA.E4   lune4   log file (ASCII)
     User:      lunti   user terminal input
        lunto   user terminal output

Files are unformatted direct-access, except that ASCII files are
formatted sequential access.
------------------------------------------------------------------------------

>>>>>>  EMF11C notes:

11/9/86  G.Miller
        Created EMF11C from EMF11B.

Compile via:  FORT/EXTEND EMF11C
Link via:     LINK EMF11C  or  LINK EMF11C,VIRP

------------------------------------------------------------------------------

>>>>>>  Main changes from EMF11B to EMF11C:

1)  A single input file is assumed, rather than a separate file for each
    trial block.  Output is to .F% files rather than .E% files.
2)  Expanded use of parameter in spec file, 'HEM'.  Setting 'HEM' to
    '-1' disables blink detection and correction.  When EMF11C is
    used to correct EEG data for HEM artifact, after EMF11B has been
    used to correct HEM and EEG for VEM artifact, this is appropriate.
3)  In the OPEN statements, the MAXREC keyword is used instead of
    the INITIALSIZE keyword for 3 of the 4 output files.  This is more
    general, since the needed output file size is computed from the
    program's parameters at run-time.

------------------------------------------------------------------------------

>>>>>>  EMF11C file usage:

     Input:     INP:zzzzzA.E3   lunin   input data
        DK:EMF11C.SPC   lunsp   specifications (ASCII)

     Output:    OUT:zzzzzA.F1   lune1   corrected averages across All blocks
        OUT:zzzzzA.F2   lune2   uncorrected averages
        OUT:zzzzzA.F3   lune3   corrected single trials
        OUT:zzzzzA.F4   lune4   log file (ASCII)

     User:      lunti   user terminal input
        lunto   user terminal output

Files are unformatted direct-access, except that ASCII files are
formatted sequential access.

------------------------------------------------------------------------------

>>>>>>  EMF11C parameter usage:

Parameters (first 10 read from specification file on lunsp):

            name                        content
            -----------------------------------------------------
            NPTS# of data points per channel per trial

            NIDS# of IDs per channel per trial

            NCHAN       # of channels total

            NCHAN1      # of EEG channels
        Note that program temporarily increments
        NCHAN1 to include EEG + EOG in program

            HEM type & location of EOG channel(s); type of correction
        -1=disable blink detection & correction
           (e.g., this is the second run of EMF11C
           on a data set -- the first run removed VEM;
           the second will remove HEM, which does not
           involve blink detection & removal)
         0=no HEM in data set
         1=HEM is the second channel; exclude the HEM
           channel from flat-spot checking in
           subroutine BADTRL
         2=HEM is the first channel; swap it with the
           second channel in subroutine STIN so that
           VEM is first & HEM is second; exclude the
           HEM channel from flat-spot checking in
           subroutine BADTRL

            MSECA/D digitizing period in msec

            EOGSEN      EOG digitization sensitivity
        (if 0, EOGSEN obtained from IDs*
        if <0, EOGSEN obtained from 1/IDs*
        if >0, used as specified)
* To permit data to remain in INTEGER*2 format but sensitivity factor to be
  in REAL*4 format, two INTEGER words are used to store the factor, if EOGSEN
  is .LE. 0.  The last 2 IDs are used -- ID(NIDS-1) and ID(NIDS).

            OUTOPT      single trial option
        0=use trials selected into a non-garbage
          bin and output all such trials
        1=use trials selected into a non-garbage
          bin; don't output corrected single trials
        2=use all good trials; unselected trials
          are put into bin # 'maxbin'; output all
          corrected single trials

            IDALL       extended IDs option (1=no, 0=yes, else=no)

            ADCRIT      off-scale A/D value
        Parameter LAB replaced by parameter ADCRIT
        (If 0, ADCRIT obtained from ID(NIDS-3))

            NBIN# of bins found
        NBIN is # experimental conditions -- event-
        related activity is computed and subtracted
        separately from each bin; correction factors
        are computed by combining all residual data
        across bins.  NBIN is controlled by the .SPC
        file and is determined by program in
        subroutine AVOUT (it's used only there and
        in SELECT).

            MAXBIN      max # bins handled (including 1 for garbage)

            MAXTRL      max # of trials handled (NOUT size)

------------------------------------------------------------------------------

>>>>>>  EMF11C specification file example (see SELECT and CHECK subroutines
for explanation of bin selection logic re: the last 8 lines in this example):

       100      NPTS    # data points following IDs     EMF11C.SPC
        30      NIDS    # IDs   11/2/86
         6      NCHAN   # A/D channels (EOG first)
         5      NCHAN1  # EEG to correct (right after EOG)
        -1      HEM     HEM channel -1,0,1,2=(see documentation)
         8      MSEC    # msec/sample (A/D sample period)
      -1.0000000EOGSEN  EOG sensitivity -1,0,>0=(see documentation)
         0      OUTOPT  0=write single trials 1=don't 2=(see docum.)
         0      IDALL   0=keep IDs 1=don't
    2047.0000000ADCRIT  off-scale criterion, 0=use ID(NIDS-3)
  1 21 1.0      bin 1:  ID21, stimtype 1 (low pitch)
  1 24 1.0      bin 1:  ID24, count
  2 21 2.0      bin 2:  ID21, stimtype 2 (high pitch)
  2 24 1.0      bin 2:  ID24, count
  3 21 1.0      bin 3:  ID21, stimtype 1 (low pitch)
  3 24 2.0      bin 3:  ID24, ignore
  4 21 2.0      bin 4:  ID21, stimtype 2 (high pitch)
  4 24 2.0      bin 4:  ID24, ignore
------------------------------------------------------------------------------

>>>>>>  EMF11D notes:

4/4/87  G.Miller
        Created EMF11D from 11/9/86 version of EMF11C.

5/8-7/17/87 C.Yee & G.Miller
        Minor changes and debugging.

Compile via:  FORT/EXTEND EMF11D
Link via:     LINK EMF11D  or  LINK EMF11D,VIRP

------------------------------------------------------------------------------

>>>>>>  Main changes from EMF11C to EMF11D:

1) Customized for experiment DT1.
2) Have user specify more of filename template (filenames no longer forced to
   "T1xnn").
3) Changed a few variable names back to the more generic names used in EMF11B.
4) Add coding of trial type into ID(18) within subroutine 'stin' -- for
   auditory trials, 1; for visual trials, 2.

------------------------------------------------------------------------------

>>>>>>  EMF11D file usage:

     Input:     INP:zzzzzA.in   lunin   input data
        DK:EMF11D.SPC   lunsp   specifications (ASCII)

     Output:    OUT:zzzzzA.o1   lune1   corrected averages across All blocks
        OUT:zzzzzA.o2   lune2   uncorrected averages
        OUT:zzzzzA.o3   lune3   corrected single trials
        OUT:zzzzzA.o4   lune4   log file (ASCII)

     User:      lunti   user terminal input
        lunto   user terminal output

User enters strings "zzzzz", "in", and "o" at run-time.
Files are unformatted direct-access, except that ASCII files are
formatted sequential access.

------------------------------------------------------------------------------

>>>>>>  EMF11D parameter usage:

Identical to EMF11C.

------------------------------------------------------------------------------

>>>>>>  EMF11D specification file example

       125      NPTS    # data points following IDs     EMF11D.SPC
        20      NIDS    # IDs   6/14/87
         4      NCHAN   # A/D channels (EOG first)
         3      NCHAN1  # EEG to correct (right after EOG)
         0      HEM     HEM channel -1,0,1,2=(see documentation)
         8      MSEC    # msec/sample (A/D sample period)
       1.5      EOGSEN  EOG sensitivity 0=use ID(NIDS-4) -1=compute
         0      OUTOPT  0=write single trials 1=don't 2=(see docum.)
         0      IDALL   0=keep IDs 1=don't
         0      ADCRIT  off-scale criterion, 0=use ID(NIDS-3)
  1  9 1.0      bin 1:   ID9,  TBtype 1  (Single - Aud)
  1 16 1.0      bin 1:   ID16, correct    (score = 1)
  2  9 2.0      bin 2:   ID9,  TBtype 2  (Single - Vis/1)
  2 16 1.0      bin 2:   ID16, corr detec (score = 1)
  3  9 2.0      bin 3:   ID9,  TBtype 2  (Single - Vis/1)
  3 16 2.0      bin 3:   ID16, corr rejec (score = 2)
  4  9 3.0      bin 4:   ID9,  TBtype 3  (Single - Vis/2)
  4 16 1.0      bin 4:   ID16, corr detec (score = 1)
  5  9 3.0      bin 5:   ID9,  TBtype 3  (Single - Vis/2)
  5 16 2.0      bin 5:   ID16, corr rejec (score = 2)
  6  9 4.0      bin 6:   ID9,  TBtype 4  (Dual - Vis/1)
  6 16 1.0      bin 6:   ID16, corr detec (score = 1)
  6 18 2.0      bin 6:   ID18, visual trial
  7  9 4.0      bin 7:   ID9,  TBtype 4  (Dual - Vis/1)
  7 16 2.0      bin 7:   ID16, corr rejec (score = 2)
  7 18 2.0      bin 7:   ID18, visual trial
  8  9 4.0      bin 8:   ID9,  TBtype 4  (Dual - Vis/1)
  8 16 1.0      bin 8:   ID16, correct    (score = 1)
  8 18 1.0      bin 8:   ID18, auditory trial
  9  9 5.0      bin 9:   ID9,  TBtype 5  (Dual - Vis/2)
  9 16 1.0      bin 9:   ID16, corr detec (score = 1)
  9 18 2.0      bin 9:   ID18, visual trial
 10  9 5.0      bin 10:  ID9,  TBtype 5  (Dual - Vis/2)
 10 16 2.0      bin 10:  ID16, corr rejec (score = 2)
 10 18 2.0      bin 10:  ID18, visual trial
 11  9 5.0      bin 11:  ID9,  TBtype 5  (Dual - Vis/2)
 11 16 1.0      bin 11:  ID16, correct    (score = 1)
 11 18 1.0      bin 11:  ID18, auditory trial
------------------------------------------------------------------------------

>>>>>>  EMF11P notes:

10/4/87  G.Miller
        Created EMF11P from 7/17/87 version of EMF11D.
1/24/88  G. Miller
        Polishing for public release.
2/28/88  G. Miller
        Final preparation for public release.

Compile via:  FORT/EXTEND EMF11P
Link via:     LINK EMF11P  or  LINK EMF11P,VIRP

------------------------------------------------------------------------------

>>>>>>  Main changes from EMF11D to EMF11P:

1. All files are now formatted sequential (ASCII), to facilitate use on other
   computers.
2. Reduced array 'NOUT' and .o3 output file 'OPEN' statement to handle 25
   instead of 1500 trials.
------------------------------------------------------------------------------

>>>>>>  EMF11P file usage:

     Input:     INP:zzzzzA.in   lunin   input data
        DK:EMF11P.SPC   lunsp   specifications (ASCII)

     Output:    OUT:zzzzzA.o1   lune1   corrected averages across All blocks
        OUT:zzzzzA.o2   lune2   uncorrected averages
        OUT:zzzzzA.o3   lune3   corrected single trials
        OUT:zzzzzA.o4   lune4   log file (ASCII)

     User:      lunti   user terminal input
        lunto   user terminal output

User enters strings "zzzzz", "in", and "o" at run-time.
Files are ASCII = formatted sequential access.

------------------------------------------------------------------------------

>>>>>>  EMF11P parameter usage:

Identical to EMF11C and EMF11D.

------------------------------------------------------------------------------

>>>>>>  EMF11P specification file example

       125      NPTS    # data points following IDs     EMF11P.SPC
        20      NIDS    # IDs   2/27/88
         4      NCHAN   # A/D channels (EOG first)
         3      NCHAN1  # EEG to correct (right after EOG)
         0      HEM     HEM channel -1,0,1,2=(see documentation)
         8      MSEC    # msec/sample (A/D sample period)
       1.5      EOGSEN  EOG sensitivity 0=use ID(NIDS-4) -1=compute
         0      OUTOPT  0=write single trials 1=don't 2=(see docum.)
         0      IDALL   0=keep IDs 1=don't
         0      ADCRIT  off-scale criterion, 0=use ID(NIDS-3)
  1  9  1.0     bin 1:   ID9,  TBtype 1  (Single Task - Aud)
  1 13 11.0     bin 1:   ID13, Stimulus Code '11'
  1 16  1.0     bin 1:   ID16, Correct    (score = 1)
  2  9  1.0     bin 2:   ID9,  TBtype 1  (Single Task - Aud)
  2 13 11.0     bin 2:   ID13, Stimulus Code '11'
  2 16  2.0     bin 2:   ID16, Incorrect  (score = 2)
  3  9  1.0     bin 3:   ID9,  TBtype 1  (Single Task - Aud)
  3 13 12.0     bin 3:   ID13, Stimulus Code '12'
  3 16  1.0     bin 3:   ID16, Correct    (score = 1)
  4  9  1.0     bin 4:   ID9,  TBtype 1  (Single Task - Aud)
  4 13 12.0     bin 4:   ID13, Stimulus Code '12'
  4 16  2.0     bin 4:   ID16, Incorrect  (score = 2)

