diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..2d407a9
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+/output
+/data
+/ignore
+*.DS_Store
\ No newline at end of file
diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100644
index 0000000..e72bfdd
--- /dev/null
+++ b/LICENSE.txt
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc.
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+
+ Copyright (C)
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ Copyright (C)
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+.
\ No newline at end of file
diff --git a/README.txt b/README.txt
new file mode 100644
index 0000000..a2f4765
--- /dev/null
+++ b/README.txt
@@ -0,0 +1,29 @@
+If this code is used in a publication, please cite the manuscript:
+"CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, GA Worrell, KJ Miller, and D Hermes.
+A preprint is available currently at doi: ___
+
+Correspondence:
+H Huang: huang.harvey@mayo.edu; D Hermes: hermes.dora@mayo.edu
+
+*****
+
+DEPENDENCIES
+
+- mnl_ieegBasics: https://github.com/MultimodalNeuroimagingLab/mnl_ieegBasics
+- matmef: https://github.com/MultimodalNeuroimagingLab/matmef
+- mnl_seegview: https://github.com/MultimodalNeuroimagingLab/mnl_seegview
+- SPM12 (as dependency for mnl_seegview): https://www.fil.ion.ucl.ac.uk/spm/software/spm12/
+- vistasoft: https://github.com/vistalab/vistasoft
+- Freesurfer v7: https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall
+- GIfTI: https://github.com/gllmflndn/gifti
+
+*****
+
+USAGE
+
+a. The data used in this analysis will be available upon manuscript publication on OpenNeuro, in BIDS format. Please download all the data and copy into a folder named "data", in the root level code directory. In the meantime before data release, all results pertaining to simulated CCEPs can still be reproduced by the code alone.
+
+b. Pial, cortical, and subcortical segmentations for each subject were obtained using Freesurfer v7. These were used to localize each stimulation site to a particular tissue type. The relevant Freesurfer outputs for each subject are located in the data subdirectory: data/derivatives/freesurfer
+
+c. All custom analyses were performed in MATLAB R2023a. Step-by-step code blocks and instructions to generate all manuscript figures and results are in "main.m". Please open "main.m" and follow the code sections contained therein.
\ No newline at end of file
diff --git a/applyCARLARealCCEPs.m b/applyCARLARealCCEPs.m
new file mode 100644
index 0000000..e179bfb
--- /dev/null
+++ b/applyCARLARealCCEPs.m
@@ -0,0 +1,217 @@
+%% This script generates outputs for CARLA applied to single stim sites in real data (Fig 5B-F, Fig 6C-F)
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+%% Load all data for selected subject
+
+% 'sub' and 'run' variables should have been defined in main script
+
+elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub));
+elecs = ieeg_readtableRmHyphens(elecsPath);
+
+mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run));
+channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run));
+eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run));
+annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run));
+
+% Load CCEP data
+mef = ccep_PreprocessMef(mefPath, channelsPath, eventsPath);
+mef.loadMefAll;
+mef.highpass;
+mef.loadMefTrials([-1, 1]);
+
+mef.plotOutputs(1, [], 200);
+
+tt = mef.tt;
+srate = mef.srate;
+channels = mef.channels;
+events = mef.evts;
+data = mef.data;
+
+sites = groupby(events.electrical_stimulation_site);
+
+outdir = fullfile('output', 'realCCEPs');
+mkdir(outdir);
+
+try
+ annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t');
+ annots = annots.status_description;
+ assert(length(annots) == height(events), 'Error: length of event annotations does not match events table height');
+catch
+ warning('Unable to load event annotation for events');
+end
+
+% show the seizure onset zones which are excluded from stimulation
+SOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ'));
+disp('SOZs:');
+disp(SOZ');
+
+clear *Path*
+
+%% Extract data from site of interest, NaN out individual bad trials
+
+% 'site' variable should have been defined in main script
+
+idxesSite = sites{strcmp(site, sites(:, 1)), 2}; % indices corresponding to this stim site
+nTrials = length(idxesSite);
+
+goodChs = strcmp(channels.type, 'SEEG') & strcmp(channels.status, 'good') & ~ismember(upper(channels.name), getNeighborChs(split(site,'-'), 2)); % exclude sites 2 away from current
+
+dataSite = data(goodChs, :, idxesSite);
+chsSite = channels.name(goodChs);
+
+% remove bad trials and set channels to nan in individual trials, according to annotations
+annotsSite = annots(idxesSite); % get annotations corresponding to this stim site
+dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite);
+
+figure('Position', [200, 200, 1200, 1200]);
+subplot(1, 2, 1); % first half of channels
+ieeg_plotTrials(tt, mean(dataSite(1:floor(length(chsSite)/2), :, :), 3, 'omitnan')', 200, chsSite(1:floor(length(chsSite)/2)));
+xlim([-0.5, 1]);
+subplot(1, 2, 2); % first half of channels
+ieeg_plotTrials(tt, mean(dataSite(ceil(length(chsSite)/2):end, :, :), 3, 'omitnan')', 200, chsSite(ceil(length(chsSite)/2):end));
+xlim([-0.5, 1]);
+
+%% Apply CARLA and plot variance, zminmean plots
+
+rng('default');
+
+nanChs = any(isnan(dataSite), [2, 3]); % any channels with nan in any trial
+
+dataSiteNoNan = dataSite(~nanChs, :, :);
+chsSiteNoNan = chsSite(~nanChs);
+
+[Vcar, CAR, stats] = CARLA(tt, dataSiteNoNan, srate, true);
+nCAR = length(stats.chsUsed);
+
+cmSens = [1, 165/255, 0]; % orange color for max
+
+% ZminMean plot
+figure('Position', [200, 200, 600, 300]); hold on
+errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k.-', 'MarkerSize', 10, 'CapSize', 1); % SD of bootstrapped means
+plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens);
+xlim([0, length(chsSiteNoNan)+1]); ylim([-1.5, 0]);
+yline(0, 'Color', 'k');
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_zMinMean', sub, site)), 'svg');
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_zMinMean', sub, site)), 'png');
+
+% Variance plot
+vars = stats.vars(stats.order); % sort in order
+figure('Position', [200, 200, 600, 300]); hold on
+plot(vars, 'k.-', 'MarkerSize', 10);
+xline(nCAR + 0.5, 'Color', cmSens);
+xlim([0, length(chsSiteNoNan)+1]);
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_covar', sub, site)), 'png');
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_covar', sub, site)), 'svg');
+
+%% Plot channels sorted by CARLA order and CAR
+
+% Plot (pre-CAR) data with line noise removed, in order of variance
+dataSite2Plot = zeros(size(dataSiteNoNan));
+for tr = 1:size(dataSiteNoNan, 3)
+ dataSite2Plot(:, :, tr) = ieeg_notch(dataSiteNoNan(:, :, tr)', srate, 60)';
+end
+dataSite2Plot = dataSite2Plot(stats.order, :, :);
+
+yspace = 120;
+
+% make excluded channels gray
+cm = zeros(length(chsSiteNoNan), 3);
+cm(nCAR+1:end, :) = repmat([0.5, 0.5, 0.5], length(chsSiteNoNan)-nCAR, 1);
+
+figure('Position', [200, 200, 1200, 1200]);
+nHalf = floor(length(chsSiteNoNan)/2); % number at which to break channels into 2 columns
+subplot(1, 2, 1); % first half of channels
+ieeg_plotTrials(tt, mean(dataSite2Plot(1:nHalf, :, :), 3)', yspace, 1:nHalf, cm(1:nHalf, :));
+xlim([-0.1, 0.5]);
+subplot(1, 2, 2); % second half of channels
+ieeg_plotTrials(tt, mean(dataSite2Plot(nHalf+1:end, :, :), 3)', yspace, nHalf+1:length(chsSiteNoNan), cm(nHalf+1:end, :));
+xlim([-0.1, 0.5]);
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs', sub, site)), 'png');
+print(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs', sub, site)),'-depsc2', '-r300', '-painters');
+
+% Plot the adjusted common average itself
+figure('Position', [200, 200, 600, 200]);
+hold on
+plot(tt, CAR, 'Color', [0.5, 0.5, 0.5]);
+plot(tt, mean(CAR, 2), 'k-', 'LineWidth', 1.5);
+xlim([-0.1, 0.5]); ylim([-1200, 1200]);
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_CAR', sub, site)), 'png');
+print(gcf, fullfile(outdir, sprintf('%s_%s_CAR', sub, site)), '-depsc2', '-r300', '-painters');
+
+%% Plot all channels in one column (for Figure S1)
+
+% downsample plotting signals to make less taxing on illustrator
+dataSite2PlotDs = nan(size(dataSite2Plot, 1), length(tt)/4, size(dataSite2Plot, 3));
+for ii = 1:size(dataSite2Plot, 3)
+ dataSite2PlotDs(:, :, ii) = downsample(dataSite2Plot(:, :, ii)', 4)';
+end
+ttDs = tt(1:4:end); % downsample time vector
+
+yspace = 80;
+
+figure('Position', [200, 200, 400, 1200]);
+ys = ieeg_plotTrials(ttDs, mean(dataSite2PlotDs, 3)', yspace, 1:length(chsSiteNoNan), cm);
+xlim([-0.1, 0.5]); ylim([ys(end) - 7*yspace, ys(1) + 3*yspace]);
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs_column', sub, site)), 'png');
+print(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs_column', sub, site)), '-depsc2', '-r300', '-painters');
+
+%% Plot examples of channels before and after CARLA
+
+% chs2Plot should be defined in main script
+
+for ii = 1:length(chs2Plot)
+
+ ch = chs2Plot(ii);
+
+ fprintf('Ch %s\n', chsSiteNoNan{stats.order(ch)});
+
+ sigRaw = squeeze(dataSiteNoNan(stats.order(ch), :, :)); % only high-pass, no notch
+ sigNotch = squeeze(dataSite2Plot(ch, :, :)); % only notch, no re-reference (order already sorted)
+ sigCar = squeeze(Vcar(stats.order(ch), :, :));
+
+ % plot raw signal
+ figure('Position', [200, 200, 450, 450]);
+ subplot(3, 1, 1); hold on;
+ plot(tt, sigRaw, 'Color', [0.5, 0.5, 0.5]);
+ plot(tt, mean(sigRaw, 2), 'k-', 'LineWidth', 1.5);
+ xlim([-0.1, 0.5]); ylim([-1200, 1200]);
+
+ % plot notched signal
+ subplot(3, 1, 2); hold on;
+ plot(tt, sigNotch, 'Color', [0.5, 0.5, 0.5]);
+ plot(tt, mean(sigNotch, 2), 'k-', 'LineWidth', 1.5);
+ xlim([-0.1, 0.5]); ylim([-600, 600]);
+
+ % plot CARLA signal
+ subplot(3, 1, 3); hold on;
+ plot(tt, sigCar, 'Color', [0.5, 0.5, 0.5]);
+ plot(tt, mean(sigCar, 2), 'k-', 'LineWidth', 1.5);
+ xlim([-0.1, 0.5]); ylim([-600, 600]);
+
+ saveas(gcf, fullfile(outdir, sprintf('%s_%s_exCh%d', sub, site, ch)), 'png');
+ print(gcf, fullfile(outdir, sprintf('%s_%s_exCh%d', sub, site, ch)), '-depsc2', '-r300', '-painters');
+
+end
+
diff --git a/applyCARLARealCCEPsLoop.m b/applyCARLARealCCEPsLoop.m
new file mode 100644
index 0000000..0bc346a
--- /dev/null
+++ b/applyCARLARealCCEPsLoop.m
@@ -0,0 +1,153 @@
+%% This loops through stim sites for each subject and saves CARLA results, with fewer other outputs along the way
+% Then run fig7_summarizeCARLARealCCEPs to summarize CARLA results for Figure 7
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+%% Paths, load data for each CCEP run
+
+outdir = fullfile('output', 'realCCEPsLoop');
+
+% subjects and runs corresponding to each subject to load
+subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'};
+runs = {[1, 2, 3, 5], [12, 13, 14], [2, 3, 4], [1, 2]};
+
+cmSens = [1, 165/255, 0]; % orange color used for detecting optimum
+
+for ss = 1:length(subs)
+
+ sub = subs{ss};
+ mkdir(outdir, sub);
+
+ % iterate across runs in subject
+ for rr = 1:length(runs{ss})
+
+ run = runs{ss}(rr);
+
+ mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run));
+ channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run));
+ eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run));
+ annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run));
+ elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub));
+
+ % Load CCEP data from each subject
+ mef = ccep_PreprocessMef(mefPath, channelsPath, eventsPath);
+ mef.loadMefAll;
+ mef.highpass;
+ mef.loadMefTrials([-1, 1]);
+
+ tt = mef.tt;
+ srate = mef.srate;
+ channels = mef.channels;
+ events = mef.evts;
+ data = mef.data;
+
+ sites = groupby(events.electrical_stimulation_site);
+
+ % read manual annotations on events
+ try
+ annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t');
+ annots = annots.status_description;
+ catch
+ warning('Unable to load event annotation for events');
+ end
+
+ % check that annots match events and remove bad events from annots
+ if exist('annots', 'var')
+ eventsTemp = ieeg_readtableRmHyphens(eventsPath, 'electrical_stimulation_site', 1); % to remove bad events from annotations
+ assert(length(annots) == height(eventsTemp), 'Error: length of event annotations does not match events table height');
+ annots(~strcmpi(eventsTemp.status, 'good')) = [];
+ end
+
+ elecs = ieeg_readtableRmHyphens(elecsPath);
+ SOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ'));
+ disp('SOZs:');
+ disp(SOZ');
+
+ %% Iterate across sites in current run
+
+ for kk = 1:size(sites, 1)
+
+ site = sites{kk};
+ idxesSite = sites{kk, 2};
+
+ if any(ismember(split(site, '-'), SOZ))
+ fprintf('Skipping stim site %s in SOZ\n', site);
+ continue
+ end
+ fprintf('Current site: %s\n', site);
+
+ % extract good channels at this stim site
+ goodChs = strcmp(channels.type, 'SEEG') & strcmp(channels.status, 'good') & ~ismember(upper(channels.name), getNeighborChs(split(site,'-'), 2)); % exclude sites 2 away from current
+ dataSite = data(goodChs, :, idxesSite);
+ chsSite = channels.name(goodChs);
+
+ try
+ % remove bad trials and set channels to nan in singular trials, according to annotations
+ annotsSite = annots(idxesSite); % get annotations corresponding to this stim site
+ dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite);
+ end
+
+ % minimum floor on number of trials for any stim site
+ if size(dataSite, 3) < 8
+ fprintf('Skipping stim site %s with %d trials\n', site, size(dataSite, 3));
+ continue
+ end
+
+ %% Remove channels that have any nan trials, apply CARLA
+
+ rng('default');
+
+ nanChs = any(isnan(dataSite), [2, 3]); % any channels with nan in any trial
+ dataSiteNoNan = dataSite(~nanChs, :, :);
+ chsSiteNoNan = chsSite(~nanChs);
+
+ % calculate average interchannel R-squared at each possible adjusted CAR size
+ rsqs = crossChCorr(tt, dataSiteNoNan, srate);
+ save(fullfile(outdir, sub, sprintf('%s_%s_interChRsq.mat', sub, site)), 'rsqs');
+
+ % Perform CAR
+ [Vcar, CAR, stats] = CARLA(tt, dataSiteNoNan, srate, true);
+ nCAR = length(stats.chsUsed);
+
+ % ZminMean plot
+ f = figure('Position', [200, 200, 600, 300]); hold on
+ errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k.-', 'MarkerSize', 10, 'CapSize', 1); % SD of bootstrapped means
+ plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens);
+ xlim([0, length(chsSiteNoNan)+1]); ylim([-inf, 0]);
+ yline(0, 'Color', 'k');
+ saveas(gcf, fullfile(outdir, sub, sprintf('%s_%s_zMinMean', sub, site)), 'svg');
+ saveas(gcf, fullfile(outdir, sub, sprintf('%s_%s_zMinMean', sub, site)), 'png');
+ close(f);
+
+ % save the optimum CAR size and number of total channels to text
+ fid = fopen(fullfile(outdir, sub, sprintf('%s_%s_nCAR.txt', sub, site)), 'w');
+ fprintf(fid, '%d\n%d', nCAR, length(chsSiteNoNan));
+ fclose(fid);
+
+
+ end
+
+ end
+
+end
\ No newline at end of file
diff --git a/fig1_simCCEPComponents.m b/fig1_simCCEPComponents.m
new file mode 100644
index 0000000..b72b586
--- /dev/null
+++ b/fig1_simCCEPComponents.m
@@ -0,0 +1,182 @@
+%% Figure 1. Create simulated CCEP data and their components
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+
+srate = 4800;
+tt = -0.5:1/srate:1-1/srate;
+
+nchs = 10; % number of simulated channels
+ntrs = 12; % number of trials
+chs = arrayfun(@(x) sprintf('ch%d', x), 1:nchs, 'UniformOutput', false)';
+
+outdir = fullfile('output', 'simComponents');
+mkdir(outdir)
+
+%% i) Create data, add stim artifacts
+
+% stores the data
+V0 = zeros(length(tt), nchs);
+
+Aart = 50 + rand(nchs, 1)*6; % slightly different artifact amplitudes for each channel
+artifact = sin(2*pi*600*tt)';
+artifact(tt < 0 | tt > 0.002) = 0;
+V0 = V0 + artifact*Aart';
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V0, 100);
+xlabel('Time (s)'); ylabel('Channels');
+
+%% A) Add unique input signals at a subset of channels
+
+rng(14);
+
+chsResp = 1:4;
+
+A = 100;
+% using previous version of signal generator for this illustration to preserve backwards compatibility
+sig = genRandSigOrig(tt, length(chsResp), A);
+
+V1 = V0;
+V1(:, chsResp) = V0(:, chsResp) + sig;
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V1, 100, [], [], 'LineWidth', 1.5);
+xlabel('Time (s)'); ylabel('Channels');
+
+%% B) Add a global signal across all chs and trials when present. Here amplitude Aglobal = 0, so no global signal added
+
+rng(10);
+
+Aglobal = 0;
+
+sigGlob = genRandSigOrig(tt, 1, Aglobal);
+
+V2 = V1 + sigGlob;
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V2, 100, [], [], 'LineWidth', 1.5);
+xlabel('Time (s)'); ylabel('Channels');
+
+%% C) Add noise common to all channels (line noise and some brown noise from reference)
+
+rng('default');
+
+% add trials, so now put V in ch x time x trial format
+V3 = repmat(V2', 1, 1, ntrs);
+
+% random line noise phases for each trial, trials x harmonic (60, 120, 180)
+phLN = rand(ntrs, 3)*2*pi;
+LN = zeros(length(tt), ntrs);
+for ii = 1:ntrs
+ LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3));
+end
+
+% brown noise shared across channels (from ref electrode
+BN = cumsum(0.4*randn(2*length(tt), ntrs));
+BN = ieeg_highpass(BN, srate, true);
+BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :);
+
+noiseCommon = LN + BN;
+
+V3 = V3 + shiftdim(noiseCommon, -1);
+
+figure('Position', [200, 200, 600, 400]);
+subplot(1, 2, 1); ieeg_plotTrials(tt, V3(:, :, 1)', 100, [], [], 'LineWidth', 1.5);
+title('V at one trial');
+xlabel('Time (s)'); ylabel('Channels');
+
+subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V3, 3)', 100, [], [], 'LineWidth', 1.5);
+title('V across trials');
+xlabel('Time (s)'); ylabel('Channels');
+
+
+%% D) add random brown noise across all channels
+
+rng('default');
+
+noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it
+for ii = 1:nchs
+ noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true);
+end
+noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :);
+
+V4 = V3 + noiseRand;
+
+figure('Position', [200, 200, 600, 400]);
+subplot(1, 2, 1); ieeg_plotTrials(tt, V4(:, :, 1)', 100, [], [], 'LineWidth', 1.5);
+title('V at one trial');
+xlabel('Time (s)'); ylabel('Channels');
+
+subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V4, 3)', 100, [], [], 'LineWidth', 1.5);
+title('V across trials');
+xlabel('Time (s)'); ylabel('Channels');
+
+saveas(gcf, fullfile(outdir, 'V4'), 'svg');
+saveas(gcf, fullfile(outdir, 'V4'), 'png');
+
+
+%% Save construction components for channels of interest
+
+% Plot select channels for figure (3 and 8)
+chs2Plot = [3, 8];
+
+for ii = 1:length(chs2Plot)
+
+ ch = chs2Plot(ii);
+
+ ylims = [-120, 120];
+
+ % plot full simulation of all trials at data
+ figure('Position', [200, 200, 600, 300]); hold on
+ plot(tt, squeeze(V4(ch, :, :)), 'Color', [0.5, 0.5, 0.5]);
+ plot(tt, mean(squeeze(V4(ch, :, :)), 2), 'k', 'LineWidth', 1.5);
+ hold off
+ ylim(ylims);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'svg');
+
+ % plot the true signal
+ figure('Position', [200, 200, 200, 300]);
+ plot(tt, V1(:, ch), 'k', 'LineWidth', 2);
+ ylim(ylims); xlim([0, 0.5]);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'svg');
+
+ % plot examples of the common noise
+ figure('Position', [200, 200, 200, 300]);
+ ieeg_plotTrials(tt, noiseCommon(:, 1:3), 50, [], [0.5, 0.5, 0.5]); % plot at one trial
+ xlim([0, 0.5]);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'png');
+ saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'svg');
+
+ % plot the random noise
+ figure('Position', [200, 200, 200, 300]); hold on
+ plot(tt, squeeze(noiseRand(ch, :, :)), 'Color', [0.5, 0.5, 0.5]);
+ hold off
+ xlim([0, 0.5]); ylim(ylims);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'svg');
+end
\ No newline at end of file
diff --git a/fig2_CARLAonSingleTrial.m b/fig2_CARLAonSingleTrial.m
new file mode 100644
index 0000000..5d01df3
--- /dev/null
+++ b/fig2_CARLAonSingleTrial.m
@@ -0,0 +1,144 @@
+%% Figure 2: Methods illustration of how CARLA works, on single-trial simulated CCEP data
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+%% A) Create data, add stim artifacts and responses to the last 2 channels
+
+srate = 4800;
+tt = -0.5:1/srate:1-1/srate;
+
+rng(50);
+
+nchs = 6; % number of simulated channels
+ntrs = 1; % number of trials
+chs = arrayfun(@(x) sprintf('ch%d', x), 1:nchs, 'UniformOutput', false)';
+
+% stores the data
+V0 = zeros(length(tt), nchs);
+
+Aart = 50 + rand(nchs, 1)*5; % slightly different artifact amplitudes for each channel
+artifact = sin(2*pi*600*tt)';
+artifact(tt < 0 | tt > 0.002) = 0;
+V0 = V0 + artifact*Aart';
+
+chsResp = [nchs-1, nchs];
+
+sig = genRandSig(tt, length(chsResp), 150);
+
+V1 = V0;
+V1(:, chsResp) = V0(:, chsResp) + sig;
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V1, 80, [], [], 'LineWidth', 1.5);
+xlabel('Time (s)'); ylabel('Channels');
+
+
+%% B) Add common noise (line noise and reference brown noise) to each trial
+
+rng('default');
+
+% no line noise in this demonstration, pretend that this has been removed by filtering
+% Create line noise: random noise phase for each harmonic (60, 120, 180)
+%phNoise = rand(1, 3)*2*pi;
+%LN = sin(2*pi*60*tt - phNoise(1))' + 0.3*sin(2*pi*120*tt - phNoise(2))' + 0.1*sin(2*pi*180*tt - phNoise(3))';
+
+% Create some brown noise and low-pass as common signal
+brownNoiseRef = cumsum(0.4*randn(2*length(tt), 1));
+brownNoiseRef = ieeg_highpass(brownNoiseRef, srate, true);
+brownNoiseRef = brownNoiseRef((0.5*length(tt)+1) : 1.5*length(tt));
+
+V3 = V1 + brownNoiseRef;
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V3, 80, [], [], 'LineWidth', 1.5);
+xlabel('Time (s)'); ylabel('Channels');
+
+
+%% C) add random brown noise to each channel
+
+rng(4);
+
+noiseRand = cumsum(0.4*randn(2*length(tt), nchs)); % give double the number of time points so we can highpass it
+noiseRand = ieeg_highpass(noiseRand, srate, true);
+
+noiseRand = noiseRand((0.5*length(tt)+1) : 1.5*length(tt), :);
+
+V4 = V3 + noiseRand;
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V4, 80, [], [], 'LineWidth', 1.5);
+xlim([-0.1, 0.5])
+xlabel('Time (s)'); ylabel('Channels');
+
+%% Apply CARLA and show outputs
+
+rng('default');
+
+[Vcar, CAR, stats] = CARLA(tt, V4, srate);
+
+outdir = fullfile('output', 'simulationTrial');
+mkdir(outdir);
+
+% Plot in increasing order of variance (Figure 2A)
+V4sorted = V4(:, stats.order);
+figure('Position', [200, 200, 200, 500]); ieeg_plotTrials(tt, V4sorted, 100, [], [], 'LineWidth', 1.5);
+xlabel('Time (s)'); ylabel('Channels');
+xlim([-0.1, 0.5]);
+xline(0.01, 'Color', 'r');
+xline(0.3, 'Color', 'r');
+saveas(gcf, fullfile(outdir, 'trialExample'), 'svg');
+saveas(gcf, fullfile(outdir, 'trialExample'), 'png');
+
+% Plot zmin (Figure 2C)
+figure('Position', [200, 200, 400, 300]); hold on
+plot(stats.zMinMean, 'k-o', 'LineWidth', 1, 'MarkerFaceColor', 'k');
+xlim([1, 6]); ylim([-1.5, 0]);
+saveas(gcf, fullfile(outdir, 'trialExampleRmin'), 'svg');
+saveas(gcf, fullfile(outdir, 'trialExampleRmin'), 'png');
+
+% Plot variance in increasing order (Figure 2D)
+vars = stats.vars(stats.order);
+figure('Position', [200, 200, 400, 250]); hold on
+plot(vars, 'k-o', 'LineWidth', 1, 'MarkerFaceColor', 'k');
+xlim([1, 6]);
+saveas(gcf, fullfile(outdir, 'trialExampleVar'), 'svg');
+saveas(gcf, fullfile(outdir, 'trialExampleVar'), 'png');
+
+%% Plot CAR-rereferenced signal at a few sizes
+
+% how many channels to use for common average to illustrate re-referenced signals (Figure 2B)
+ns = [2, 4, 5];
+
+for ii = 1:length(ns)
+
+ n = ns(ii);
+
+ CARcurr = mean(V4sorted(:, 1:n), 2);
+ V4Carn = V4sorted - mean(V4sorted(:, 1:n), 2);
+
+ figure('Position', [200, 200, 150, 500]); ieeg_plotTrials(tt, V4Carn, 100, [], [], 'LineWidth', 1.5);
+ xlabel('Time (s)'); ylabel('Channels');
+ xlim([0.01, 0.3]);
+ saveas(gcf, fullfile(outdir, sprintf('trialExampleCAR%dch', n)), 'svg');
+ saveas(gcf, fullfile(outdir, sprintf('trialExampleCAR%dch', n)), 'png');
+
+end
+
diff --git a/fig3_EPConstruct.m b/fig3_EPConstruct.m
new file mode 100644
index 0000000..030a8a0
--- /dev/null
+++ b/fig3_EPConstruct.m
@@ -0,0 +1,128 @@
+%% This script plots an example of enveloped sinusoids that add to form a simulated evoked potential (Figure 3)
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+%% Configure signal parameters
+
+rng(4);
+
+srate = 4800;
+tt = -0.5:1/srate:1-1/srate;
+
+A = 100;
+
+tau1 = 0.01 + rand*0.02; % time constant to dampen the faster sinusoid (s1)
+tau2 = 0.06 + rand*0.08;
+tauOn1 = 0.005; % time constant of subtracted exponential decay for both
+tauOn2 = 0.025; % time constant of subtracted exponential decay for both
+
+f1 = 8 + rand*4; % freq (hz) of first sinusoid (mean period = 100 ms)
+f2 = 1 + rand*2; % freq of second sinusoid, (mean period = 500 ms)
+
+ph1 = rand*2*pi; % phase of first sinusoid
+ph2 = rand*2*pi; % phase of second sinusoid
+
+outdir = fullfile('output', 'EPConstruct');
+mkdir(outdir);
+
+%% Construct and save full signal (3C)
+
+sig = A*( (exp(-tt/tau1) - exp(-tt/tauOn1)).*sin(2*pi*f1*tt - ph1) + (exp(-tt/tau2) - exp(-tt/tauOn2)).*sin(2*pi*f2*tt - ph2) );
+sig(tt < 0) = 0; % make causal
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, sig, 'k-');
+xlabel('Time (s)');
+xlim([-0.1, 0.5]);
+ylim([-80, 80]);
+
+saveas(gcf, fullfile(outdir, 'signal'), 'png');
+saveas(gcf, fullfile(outdir, 'signal'), 'svg');
+
+%% Plot faster part of signal (3A)
+
+% the direct/faster signal component
+env1 = (exp(-tt/tau1) - exp(-tt/tauOn1));
+sig1 = A*sin(2*pi*f1*tt - ph1);
+env1(tt < 0) = 0; % make causal
+sig1(tt < 0) = 0;
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, env1, 'Color', [0.5, 0.5, 0.5], 'LineWidth', 2);
+xline(tau1, 'Color', 'r'); xline(tauOn1, 'Color', 'r');
+xlim([-0.1, 0.5]);
+ylim([-1, 1]);
+saveas(gcf, fullfile(outdir, 'env1'), 'png');
+saveas(gcf, fullfile(outdir, 'env1'), 'svg');
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, sig1, 'k-');
+xlim([-0.1, 0.5]);
+ylim([-120, 120]);
+saveas(gcf, fullfile(outdir, 'sig1'), 'png');
+saveas(gcf, fullfile(outdir, 'sig1'), 'svg');
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, env1.*sig1, 'k-');
+xlim([-0.1, 0.5]);
+ylim([-80, 80]);
+saveas(gcf, fullfile(outdir, 'env1sig1'), 'png');
+saveas(gcf, fullfile(outdir, 'env1sig1'), 'svg');
+
+%% Plot slower part of signal (3B)
+
+env2 = (exp(-tt/tau2) - exp(-tt/tauOn2));
+sig2 = A*sin(2*pi*f2*tt - ph2);
+env2(tt < 0) = 0; % make causal
+sig2(tt < 0) = 0;
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, env2, 'Color', [0.5, 0.5, 0.5], 'LineWidth', 2);
+xlim([-0.1, 0.5]);
+ylim([-1, 1]);
+xline(tau2, 'Color', 'r'); xline(tauOn2, 'Color', 'r');
+saveas(gcf, fullfile(outdir, 'env2'), 'png');
+saveas(gcf, fullfile(outdir, 'env2'), 'svg');
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, sig2, 'k-');
+xlim([-0.1, 0.5]);
+ylim([-120, 120]);
+saveas(gcf, fullfile(outdir, 'sig2'), 'png');
+saveas(gcf, fullfile(outdir, 'sig2'), 'svg');
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, env2.*sig2, 'k-');
+xlim([-0.1, 0.5]);
+ylim([-80, 80]);
+saveas(gcf, fullfile(outdir, 'env2sig2'), 'png');
+saveas(gcf, fullfile(outdir, 'env2sig2'), 'svg');
+
+%% Save full signal (again), using the constructed components
+
+figure('Position', [200, 200, 200, 100]);
+plot(tt, env1.*sig1 + env2.*sig2, 'k-');
+xlabel('Time (s)');
+xlim([-0.1, 0.5]);
+ylim([-80, 80]);
diff --git a/fig7_summarizeCARLARealCCEPs.m b/fig7_summarizeCARLARealCCEPs.m
new file mode 100644
index 0000000..07c8095
--- /dev/null
+++ b/fig7_summarizeCARLARealCCEPs.m
@@ -0,0 +1,162 @@
+%% This loads the outputs from applyCARLARealCCEPsLoop to summarize the optimal CCEP cutoff size for stim sites
+% in each tissue type across the 4 subjects
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+
+% subjects and runs corresponding to each subject to load
+subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'};
+runs = {[1, 2, 3, 5], [12, 13, 14], [2, 3, 4], [1, 2]};
+
+for ss = 1:length(subs)
+
+ sub = subs{ss};
+ fprintf('Sub: %s\n', sub);
+
+ %% Obtain stim sites across all runs, number of trials
+
+ elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub));
+ elecs = ieeg_readtableRmHyphens(elecsPath);
+ SOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ'));
+
+ sites = cell(0, 2); % name, number of trials
+ for rr = 1:length(runs{ss})
+ run = runs{ss}(rr);
+ fprintf('%s run %d\n', sub, run);
+
+ eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run));
+ annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run));
+
+ events = ieeg_readtableRmHyphens(eventsPath, 'electrical_stimulation_site', 1);
+ eventsStatusOrig = events.status;
+ events(~strcmpi(events.status, 'good'), :) = [];
+
+ sitesRun = groupby(events.electrical_stimulation_site); % all sites in current run
+
+ % read manual annotations on events if it exists
+ try
+ annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t');
+ annots = strip(annots.status_description);
+ catch
+ warning('Unable to load event annotation for events');
+ clear annots
+ end
+
+ if exist('annots', 'var')
+ annots(~strcmpi(eventsStatusOrig, 'good')) = [];
+ assert(length(annots) == height(events), 'Error: length of event annotations does not match events table height');
+ end
+
+ % Remove stim sites in SOZ
+ isSoz = any(ismember(split(sitesRun(:, 1), '-', 2), SOZ), 2);
+ sitesRun(isSoz, :) = [];
+
+ % remove "all" annotated trials and calculate final number of trials
+ nTrials = cellfun(@length, sitesRun(:, 2));
+ if exist('annots', 'var')
+ for ii = 1:size(sitesRun, 1)
+ annotsCurr = annots(sitesRun{ii, 2});
+ nTrsRm = sum(strcmpi(annotsCurr, 'all'));
+ nTrials(ii) = nTrials(ii) - nTrsRm;
+ end
+ end
+
+ sites = [sites; sitesRun(:, 1), num2cell(nTrials)];
+
+ end
+
+ % ensure no stim site has less than 8 trials
+ sites = cell2table(sites, 'VariableNames', {'site', 'trials'});
+ sites(sites.trials < 8, :) = [];
+
+ %% Identify Destrieux labels for each stim site
+
+ % subject-specific Freesurfer derivatives folder
+ FSdir = fullfile(dataPath, 'derivatives', 'freesurfer', sub);
+ xyzs = ieeg_getPairXyzs(sites.site, elecs);
+ labs = ieeg_getLabelXyzDestrieux(xyzs, FSdir); % within radius of 3
+
+ labGroup = repmat({'n/a'}, length(labs), 1);
+ labGroup(contains(labs, {'Thalamus', 'Hippocampus', 'Amygdala'})) = {'subgray'};
+ labGroup(contains(labs, 'White_Matter')) = {'WM'};
+ labGroup(startsWith(labs, {'lh', 'rh'})) = {'cortgray'};
+
+ fprintf('Subcortical Gray: %d\nCortical Gray: %d\nWhite Matter: %d\nn/a: %d\n', ...
+ sum(strcmp(labGroup, 'subgray')), sum(strcmp(labGroup, 'cortgray')), sum(strcmp(labGroup, 'WM')), sum(strcmp(labGroup, 'n/a')));
+
+ sites.labs = labs;
+ sites.labGroup = labGroup;
+
+ %% load the nCAR files and plot distributions
+
+ rng('default');
+
+ outdir = fullfile('output', ''realCCEPsLoop', sub);
+
+ carSz = zeros(height(sites), 2);
+ for ii = 1:height(sites)
+ M = readmatrix(fullfile(outdir, sprintf('%s_%s_nCAR.txt', sub, sites.site{ii})));
+ carSz(ii, :) = M;
+ end
+
+ sites.carSz = carSz(:, 1);
+ sites.numChs = carSz(:, 2);
+ sites.pct = 100*carSz(:, 1)./carSz(:, 2);
+
+ % plot overall distribution
+ figure; histogram(sites.pct, 10);
+
+ % plot distributions by label group and hemisphere
+ pctByGroup = cell(6, 1);
+
+ % Right hemi
+ pctByGroup{1} = sites.pct(strcmp(labGroup, 'cortgray') & startsWith(lower(labs), 'r'));
+ pctByGroup{2} = sites.pct(strcmp(labGroup, 'subgray') & startsWith(lower(labs), 'r'));
+ pctByGroup{3} = sites.pct(strcmp(labGroup, 'WM') & startsWith(lower(labs), 'r'));
+
+ % Left hemi
+ pctByGroup{4} = sites.pct(strcmp(labGroup, 'cortgray') & startsWith(lower(labs), 'l'));
+ pctByGroup{5} = sites.pct(strcmp(labGroup, 'subgray') & startsWith(lower(labs), 'l'));
+ pctByGroup{6} = sites.pct(strcmp(labGroup, 'WM') & startsWith(lower(labs), 'l'));
+
+
+ figure('Position', [200, 200, 250, 300]); hold on
+ bar(1:3, cellfun(@mean, {[pctByGroup{1}; pctByGroup{4}], [pctByGroup{2}; pctByGroup{5}], [pctByGroup{3}; pctByGroup{6}]}), ...
+ 'FaceColor', [0.8, 0.8, 0.8]); % mean for stim sites in both hemispheres
+ hsR = jitterplot(pctByGroup(1:3), brighten([252, 141, 98]/255, -0.5));
+ set(hsR, 'Marker', '.', 'MarkerSize', 16);
+ hsL = jitterplot(pctByGroup(4:6), brighten([102, 194, 165]/255, -0.5));
+ set(hsL, 'Marker', '.', 'MarkerSize', 16);
+ ylim([0, 100]);
+ set(gca, 'XTick', 1:3, 'XTickLabels', {'cortgray', 'subgray', 'WM'});
+
+ [p, tbl, stats] = anova1(sites.pct, sites.labGroup, 'off');
+ fprintf('ANOVA p = %0.02e\n', p);
+
+ saveas(gcf, fullfile(outdir, 'summaryPctCAR'), 'png');
+ saveas(gcf, fullfile(outdir, 'summaryPctCAR'), 'svg');
+
+ writetable(sites, fullfile(outdir, 'summaryTbl.tsv'), 'Delimiter', '\t', 'FileType', 'text');
+
+end
diff --git a/fig8_summarizeInterChCorr.m b/fig8_summarizeInterChCorr.m
new file mode 100644
index 0000000..e09c8b5
--- /dev/null
+++ b/fig8_summarizeInterChCorr.m
@@ -0,0 +1,125 @@
+%% This script uses output from applyCARLARealCCEPsLoop to calculate cross-channel rsqs for each adjusted CAR method
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+
+outdir = fullfile('output', 'realCCEPsLoop');
+
+% subjects
+subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'};
+
+%% Load .mat files and CARLA results
+
+% sub, stim site, z-transformed interch correlation matrix, CARLA result
+sites = cell(0, 4);
+
+for ss = 1:length(subs)
+
+ sub = subs{ss};
+
+ mats = glob(fullfile(outdir, sub, sprintf('%s_*-*_interChRsq.mat', sub)));
+
+ for jj = 1:length(mats)
+
+ % note which site is being added currently
+ eles = split(mats{jj}, '_');
+ load(mats{jj});
+ M = readmatrix(fullfile(outdir, sub, sprintf('%s_%s_nCAR.txt', sub, eles{2})));
+
+ % 3rd dim of zs == number of recording channels + 1 (corresponding to no reference)
+ assert(size(rsqs, 3) == M(2)+1, 'Mismatch between number of channels in CARLA output and correlation output');
+
+ sites = [sites; {sub}, eles(2), {rsqs}, {M}];
+ end
+end
+
+assert(size(sites, 1) == 82, 'Error: expecting 82 total stim sites across the 4 subjects');
+
+%% Summarize mean correlation for no re-referencing, 50% covariance, CARLA, and full rereferencing
+
+func = (@(x) mean(x, 'all', 'omitnan'));
+
+% each row = stim site. Order of columns = no rereferencing, 50% rereferencing, full rereferencing, CARLA
+rsqsAll = nan(size(sites, 1), 5);
+for ii = 1:size(sites, 1)
+ rsqs = sites{ii, 3};
+ nChs = sites{ii, 4}(2); % number of recording channels
+ nCarla = sites{ii, 4}(1); % optimal carla size
+
+ rsqsAll(ii, 1) = func(rsqs(:, :, 1)); % mean correlation without referencing
+ rsqsAll(ii, 2) = func(rsqs(:, :, end)); % mean correlation at full CAR
+ rsqsAll(ii, 3) = func(rsqs(:, :, floor(nChs*0.25) + 1)); % mean correlation at 25%
+ rsqsAll(ii, 4) = func(rsqs(:, :, floor(nChs*0.5) + 1)); % mean correlation at 50%
+ rsqsAll(ii, 5) = func(rsqs(:, :, nCarla + 1)); % mean correlation using CARLA optimum
+end
+
+figure('Position', [200, 200, 200, 300]);
+boxplot(rsqsAll, 'Symbol', '.');
+xlim([0, 6]); ylim([0, 1]);
+saveas(gcf, fullfile(outdir, 'rsCompsAllsubs'), 'png');
+saveas(gcf, fullfile(outdir, 'rsCompsAllsubs'), 'svg');
+
+ps = nan(4, 1);
+for ii = 1:4
+ ps(ii) = signrank(rsqsAll(:, ii), rsqsAll(:, 5));
+end
+psBonf = ps*4;
+disp(psBonf);
+
+%% Plot histogram of pair-wise stim site differences between each reference condition and CARLA
+
+rsqsDiff = rsqsAll - rsqsAll(:, 5);
+figure('Position', [200, 200, 200, 300]);
+boxplot(rsqsDiff(:, 1:4), 'Symbol', '.');
+yline(0, 'Color', [0.2, 0.2, 0.2]);
+xlim([0, 5]); ylim([-0.4, 0.8]);
+saveas(gcf, fullfile(outdir, 'rsDiffAllsubs'), 'png');
+saveas(gcf, fullfile(outdir, 'rsDiffAllsubs'), 'svg');
+
+%% Plot matrices for sub1 stimsite 1
+
+sub = 'sub-1';
+site = 'RMO8-RMO9';
+
+idx = find(strcmp(sites(:, 1), sub) & strcmp(sites(:, 2), site));
+rmat = sites{idx, 3};
+nChs = sites{idx, 4}(2);
+nCAR = sites{idx, 4}(1);
+
+% no rereferencing, standard CAR, naive low variance CAR, CARLA
+labels = {'none', 'CAR', 'bottom25', 'bottom50', 'CARLA'};
+rmatSelect = rmat(:, :, [1, nChs+1, floor(nChs*0.25)+1, floor(nChs*0.5)+1, nCAR+1]);
+
+for ii = 1:5
+
+ figure;
+ imagesc(rmatSelect(:, :, ii), [0, 1]);
+ axis square;
+ colormap(hot);
+ colorbar;
+
+ saveas(gcf, fullfile(outdir, sprintf('rsMap_%s_%s-%s', labels{ii}, sub, site)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('rsMap_%s_%s-%s', labels{ii}, sub, site)), 'svg');
+
+end
diff --git a/figS1_compareCARvCARLA.m b/figS1_compareCARvCARLA.m
new file mode 100644
index 0000000..28dfcbd
--- /dev/null
+++ b/figS1_compareCARvCARLA.m
@@ -0,0 +1,55 @@
+%% Compare CARLA vs CAR re-referenced signals on same set of channels. Run this script as an addendum after running applyCARLARealCCEPs
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+
+%% Sort the CARLA output and calculate the standard CAR-referenced signals
+
+% sort the CARLA output channels by CARLA order
+VcarlaSorted = Vcar(stats.order, :, :);
+
+% calculate the standard common average and subtract from signals
+CARstandard = mean(dataSiteNoNan, 1, 'omitnan');
+VcarStandard = dataSiteNoNan - CARstandard;
+VcarStandardSorted = VcarStandard(stats.order, :, :);
+
+%% Zoom into the first 50 channels, compare side by side CAR vs CARLA
+
+yspace = 50;
+numChsToPlot = 50;
+
+figure('Position', [200, 200, 1000, 1400]);
+
+% the CARLA channels
+subplot(1, 2, 1);
+ieeg_plotTrials(tt, mean(VcarlaSorted(1:numChsToPlot, :, :), 3)', yspace, 1:numChsToPlot, [0, 0, 0]);
+xlim([-0.1, 0.5]);
+
+% the standard CAR channels
+subplot(1, 2, 2);
+ieeg_plotTrials(tt, mean(VcarStandardSorted(1:numChsToPlot, :, :), 3)', yspace, 1:numChsToPlot, [0, 0, 0]);
+xlim([-0.1, 0.5]);
+
+saveas(gcf, fullfile(outdir, sprintf('%s_%s_sortedChsCARvCARLA_1to%d', sub, site, numChsToPlot)), 'png');
+print(gcf, fullfile(outdir, sprintf('%s_%s_sortedChsCARvCARLA_1to%d', sub, site, numChsToPlot)),'-depsc2', '-r300', '-painters');
diff --git a/figS3_simCCEPComponentsGlobalSig.m b/figS3_simCCEPComponentsGlobalSig.m
new file mode 100644
index 0000000..69ce024
--- /dev/null
+++ b/figS3_simCCEPComponentsGlobalSig.m
@@ -0,0 +1,191 @@
+%% Figure S2. Create simulated CCEP data and their components, like Figure 1, but this time containing a global signal
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+
+srate = 4800;
+tt = -0.5:1/srate:1-1/srate;
+
+nchs = 10; % number of simulated channels
+ntrs = 12; % number of trials
+chs = arrayfun(@(x) sprintf('ch%d', x), 1:nchs, 'UniformOutput', false)';
+
+outdir = fullfile('output', 'simComponentsGlobal');
+mkdir(outdir)
+
+%% i) Create data, add stim artifacts
+
+% stores the data
+V0 = zeros(length(tt), nchs);
+
+Aart = 50 + rand(nchs, 1)*6; % slightly different artifact amplitudes for each channel
+artifact = sin(2*pi*600*tt)';
+artifact(tt < 0 | tt > 0.002) = 0;
+V0 = V0 + artifact*Aart';
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V0, 100);
+xlabel('Time (s)'); ylabel('Channels');
+
+%% A) Add unique input signals at a subset of channels
+
+rng(14);
+
+chsResp = 1:4;
+
+A = 100;
+% using previous version of signal generator for this illustration to preserve backwards compatibility
+sig = genRandSigOrig(tt, length(chsResp), A);
+
+V1 = V0;
+V1(:, chsResp) = V0(:, chsResp) + sig;
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V1, 100, [], [], 'LineWidth', 1.5);
+xlabel('Time (s)'); ylabel('Channels');
+
+%% B) Add a global noise across all chs and trials, e.g. from reference.
+
+rng(10);
+
+Aglobal = 20; % same amplitude of noise on all channels
+
+sigGlob = genRandSigOrig(tt, 1, Aglobal);
+
+V2 = V1 + sigGlob;
+
+figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V2, 100, [], [], 'LineWidth', 1.5);
+xlabel('Time (s)'); ylabel('Channels');
+
+%% C) Add noise common to all channels (line noise and some brown noise from reference)
+
+rng('default');
+
+% add trials, so now put V in ch x time x trial format
+V3 = repmat(V2', 1, 1, ntrs);
+
+% random line noise phases for each trial, trials x harmonic (60, 120, 180)
+phLN = rand(ntrs, 3)*2*pi;
+LN = zeros(length(tt), ntrs);
+for ii = 1:ntrs
+ LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3));
+end
+
+% brown noise shared across channels (from ref electrode
+BN = cumsum(0.4*randn(2*length(tt), ntrs));
+BN = ieeg_highpass(BN, srate, true);
+BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :);
+
+noiseCommon = LN + BN;
+
+V3 = V3 + shiftdim(noiseCommon, -1);
+
+figure('Position', [200, 200, 600, 400]);
+subplot(1, 2, 1); ieeg_plotTrials(tt, V3(:, :, 1)', 100, [], [], 'LineWidth', 1.5);
+title('V at one trial');
+xlabel('Time (s)'); ylabel('Channels');
+
+subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V3, 3)', 100, [], [], 'LineWidth', 1.5);
+title('V across trials');
+xlabel('Time (s)'); ylabel('Channels');
+
+
+%% D) add random brown noise across all channels
+
+rng('default');
+
+noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it
+for ii = 1:nchs
+ noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true);
+end
+noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :);
+
+V4 = V3 + noiseRand;
+
+figure('Position', [200, 200, 600, 400]);
+subplot(1, 2, 1); ieeg_plotTrials(tt, V4(:, :, 1)', 100, [], [], 'LineWidth', 1.5);
+title('V at one trial');
+xlabel('Time (s)'); ylabel('Channels');
+
+subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V4, 3)', 100, [], [], 'LineWidth', 1.5);
+title('V across trials');
+xlabel('Time (s)'); ylabel('Channels');
+
+saveas(gcf, fullfile(outdir, 'V4'), 'svg');
+saveas(gcf, fullfile(outdir, 'V4'), 'png');
+
+
+%% Save construction components for channels of interest
+
+ylims = [-120, 120];
+
+% Plot select channels for figure (3 and 8)
+chs2Plot = [3, 8];
+
+for ii = 1:length(chs2Plot)
+
+ ch = chs2Plot(ii);
+
+ % plot full simulation of all trials at data
+ figure('Position', [200, 200, 600, 300]); hold on
+ plot(tt, squeeze(V4(ch, :, :)), 'Color', [0.5, 0.5, 0.5]);
+ plot(tt, mean(squeeze(V4(ch, :, :)), 2), 'k', 'LineWidth', 1.5);
+ hold off
+ ylim(ylims);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'svg');
+
+ % plot the true signal
+ figure('Position', [200, 200, 200, 300]);
+ plot(tt, V1(:, ch), 'k', 'LineWidth', 2);
+ ylim(ylims); xlim([0, 0.5]);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'svg');
+
+ % plot the global noise
+ figure('Position', [200, 200, 200, 300]);
+ plot(tt, sigGlob, 'k', 'LineWidth', 2);
+ ylim(ylims); xlim([0, 0.5]);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, 'globalSig'), 'png');
+ saveas(gcf, fullfile(outdir, 'globalNoise'), 'svg');
+
+ % plot examples of the common noise
+ figure('Position', [200, 200, 200, 300]);
+ ieeg_plotTrials(tt, noiseCommon(:, 1:3), 50, [], [0.5, 0.5, 0.5]); % plot at one trial
+ xlim([0, 0.5]);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'png');
+ saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'svg');
+
+ % plot the random noise
+ figure('Position', [200, 200, 200, 300]); hold on
+ plot(tt, squeeze(noiseRand(ch, :, :)), 'Color', [0.5, 0.5, 0.5]);
+ hold off
+ xlim([0, 0.5]); ylim(ylims);
+ xlabel('Time (s)'); ylabel('Voltage (\muV)');
+ saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'svg');
+
+end
diff --git a/functions/CARLA.m b/functions/CARLA.m
new file mode 100644
index 0000000..6fd36b2
--- /dev/null
+++ b/functions/CARLA.m
@@ -0,0 +1,225 @@
+function [Vout, CAR, stats] = CARLA(tt, V, srate, sens, nboot)
+%
+% This function performs Common Average Re-referencing by Least Anticorrelation (CARLA).
+% Channels are ranked in order of increasing covariance across trials (or variance if single trial).
+% Next, channels are iteratively added into a CAR, and the mean correlation between each unre-referenced channel and all re-referenced channels is
+% calculated, for each bootstrapped mean signal. zmin denotes the mean correlation belonging, at each bootstrap, to the most anticorrelated channel.
+% The optimum CAR size corresponds to when there is least anticorrelation (when zmin takes on its least negative value)
+%
+% Vout = CARLA(tt, V); [NOT RECOMMENDED]
+% Vout = CARLA(tt, V, srate);
+% [Vout, stats] = CARLA(tt, V, srate, nboot);
+% tt = 1 x t num. Time points matching V, in seconds
+% V = n x t x k num, or n x t num. Signal data to be re-referenced, with n channels by t time points by k trials (2 dimensional if k=1).
+% srate = num. Sampling frequency of data. If not given, this will be estimated from the time points (not recommended because of potential imprecision).
+% sens = (optional) boolean. Determines whether the more sensitive cutoff will be used, corresponding to before the first significant decrease in zMinMean.
+% Ignored if V is n x t (needs trials to calculate significance). Default = false (returns global maximum in zMinMean)
+% nboot = (optional) integer. Number of bootstrapped mean signal samples to generate when calculating correlations. Default = 100. Ignored if data is n x t
+%
+% RETURNS:
+% Vout = n x t x k num, or n x t num. Re-referenced signal data.
+% CAR = t x k num, or 1 x t num. The CAR at each trial
+% stats = struct, containing fields:
+% chsUsed = m x 1 num, m < n. The sorted n channel numbers used in the optimum CAR
+% vars = n x 1 num. Trial covariance of each channel (variance if K = 1)
+% order = n x 1 num. The channels sorted in order of increasing covariance (variance if K = 1)
+% zMin = n x n x nboot num, or n x n num if k = 1. The z (Fisher z-transformed Pearson's r) for the channel with the most negative
+% average z, at each bootstrap. Dimensions are (z to each channel) x (number of channels in CAR subset) x (bootstrap)
+% zMinMean = n x nboot num, or 1 x n num if k = 1. Mean zMin across the first dimension (averaged across all "target" rereferenced channels)
+% for each bootstrap if V was 3 dimensional
+%
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+% 2023/04/27 by Harvey Huang
+%
+ % Ensure ch x time points structure if k = 1
+ if ismatrix(V)
+ if size(V, 1) == length(tt), V = V'; end
+ end
+
+ if nargin < 5 || isempty(nboot), nboot = 100; end
+ if nargin < 4 || isempty(sens), sens = false; end
+
+ % Estimate sampling frequency from time points if necessary (not recommended)
+ if nargin < 3 || isempty(srate)
+ srate = (length(tt) - 1) / (max(tt) - min(tt));
+ warning('sampling frequency was estimated to be %0.02f', srate);
+ end
+
+ nChs = size(V, 1);
+ nTrs = size(V, 3);
+ assert(size(V, 2) == length(tt), 'Error: second dimension does not match tt. Data should be channels x timepoints x trials');
+
+ stats = struct();
+
+ % Notch filter signals to remove line noise. This is not used in the final output
+ Vclean = V;
+ for ff = 60:60:180
+ dNotch = designfilt('bandstopiir', 'FilterOrder', 4, ...
+ 'DesignMethod', 'butter', ...
+ 'HalfPowerFrequency1', ff-2, ...
+ 'HalfPowerFrequency2', ff+2, ...
+ 'SampleRate', srate);
+ for ii = 1:nTrs % filter per trial
+ Vclean(:, :, ii) = filtfilt(dNotch, Vclean(:, :, ii)')';
+ end
+ end
+
+ % extract data on response interval to perform analysis on
+ twin = [0.01, 0.3]; % should this be an input parameter? For now keep hardcoded
+ Vseg = Vclean(:, tt >= twin(1) & tt <= twin(2), :);
+
+ % Rank channels. If single trial rank based on increasing variance. If more than 1 trial, rank by covariance (self-self excluded). Better estimator of consistent structure
+ if nTrs == 1
+ stats.vars = var(Vseg, 0, 2);
+ else
+ stats.vars = nan(nChs, 1);
+ for ii = 1:nChs
+ covCurr = cov(squeeze(Vseg(ii, :, :)));
+ stats.vars(ii) = mean(covCurr(logical(triu(ones(size(covCurr)), 1))), 'all');
+ end
+ end
+ [~, stats.order] = sort(stats.vars);
+
+ % channels x CARsize x bootstrapped sample
+ if nTrs == 1
+ stats.zMin = nan(nChs, nChs); % no bootstrapping if only single trial
+ else
+ stats.zMin = nan(nChs, nChs, nboot);
+ end
+
+ % Pull subset U (order(1:ii)) from V, growing U each time by 1, to construct CAR
+ for ii = 2:nChs % min CAR size is 2
+
+ % Vseg after subtracting the putative CAR on subset U
+ VsegReref = Vseg - mean(Vseg(stats.order(1:ii), :, :), 1);
+
+ if nTrs == 1 % no bootstrapping, test correlation at the single trial level
+ Useg = Vseg(stats.order(1:ii), :)'; % UNREREFERENCED channels in subset to correlate with the rereferenced signals
+ UsegReref = VsegReref(stats.order(1:ii), :)'; % rereferenced signals in subset
+
+ r = corr(Useg, UsegReref); % rows = which unrereferenced channel was used.
+ r(1:(ii+1):end) = nan; % omit self-self correlations
+ z = atanh(r); % fisher z-transform
+
+ % choose "most responsive" channel with greatest anticorrelation, on average
+ [~, kkMost] = min(mean(z, 2, 'omitnan'));
+ stats.zMin(stats.order(1:ii), ii) = z(kkMost, :);
+
+ continue;
+
+ end
+
+ % if nTrs > 1
+ for jj = 1:nboot
+
+ % bootstrap to calculate mean signal
+ inds = datasample(1:nTrs, nTrs);
+ Useg = mean(Vseg(stats.order(1:ii), :, inds), 3)';
+ UsegReref = mean(VsegReref(stats.order(1:ii), :, inds), 3)';
+
+ r = corr(Useg, UsegReref); % rows = which unrereferenced signal was used.
+ r(1:(ii+1):end) = nan; % omit self-self correlations
+ z = atanh(r); % fisher z-transform
+
+ % choose "most responsive" channel with greatest anticorrelation, on average, across channels, for this bootstrapped sample
+ [~, kkMost] = min(mean(z, 2, 'omitnan'));
+ stats.zMin(stats.order(1:ii), ii, jj) = z(kkMost, :);
+
+ end
+
+ end
+
+ % calculate mean across the unreferenced target channels
+ stats.zMinMean = mean(stats.zMin, 1, 'omitnan');
+
+ if sens && nTrs > 1 % find sensitive optimum
+
+ nMin = ceil(0.1*nChs); % minimum number of channels to start at, for stability. Hardcoded at 10 %
+ zMMxTrs = mean(stats.zMinMean(1, :, :), 3); % mean ZMinMean across bootstrapped samples for each n
+
+ ii = nMin;
+ while ii <= nChs
+
+ if ii == nChs % at the end. no more comparisons needed
+ nOptimum = ii;
+ break
+ end
+
+ % loop until next local maximum found
+ if zMMxTrs(ii + 1) > zMMxTrs(ii)
+ ii = ii + 1;
+ continue;
+ end
+
+ zMMxTrs(1:ii-1) = nan; % set all earlier values to nan, to avoid them being found as nextGreater or as trough
+
+ % find index of the next greater zMinmean (ignores earlier nans), if it exists
+ nextGreater = find(zMMxTrs > zMMxTrs(ii), 1, 'First');
+ if isempty(nextGreater) % no more samples greater than current peak. Optimum is identical to global optimum. Done
+ nOptimum = ii;
+ break
+ end
+ zMMCurr = squeeze(stats.zMinMean(1, ii, :)); % all bootstrapped samples at current max
+ [~, idx] = min(zMMxTrs(1:nextGreater-1)); % locate Zminmean trough before the next max (ignores earlier nans)
+ zMMTrough = squeeze(stats.zMinMean(1, idx, :)); % all bootstrapped samples at trough before next greater value
+
+ % all pairwise differences in the mean correlation estimate between the trough and current peak.
+ % Note that zMinMean is the estimated correlation of the mean signal. We are testing the difference of zMinMean directly,
+ % Not the mean of zMinMean, which would be redundant and inappropriate because the bootstrapped samples are not independent.
+ [X, Y] = meshgrid(zMMTrough, zMMCurr);
+ diffs = X - Y; diffs = diffs(:);
+ conf = [-inf, prctile(diffs, 95)]; % left-tailed 95% confidence interval (H_a = difference is less than 0)
+
+ % found significant decrease. Current max is optimal.
+ if conf(2) < 0
+ nOptimum = ii;
+ break
+ end
+
+ % Difference was not significant. Start at next greater sample
+ ii = nextGreater;
+
+ end
+
+ if nOptimum == nMin, warning('Optimum CAR was detected at minimal 10% floor'); end
+
+ else
+ % Find optimum n at max (least negative) zminMean, averaged across trials.
+ [~, nOptimum] = max(mean(stats.zMinMean, 3));
+ end
+
+ % delete first single dimension if 3D, after finding nOptimum
+ stats.zMinMean = squeeze(stats.zMinMean);
+
+ % which channels were used to make CAR
+ stats.chsUsed = sort(stats.order(1:nOptimum));
+
+ CAR = mean(V(stats.chsUsed, :, :));
+
+ % Using the original (unfiltered signals)
+ Vout = V - CAR;
+
+ CAR = squeeze(CAR); % squeeze out the first dimension to make t x k, if 3 dimensional
+
+end
diff --git a/functions/antifind.m b/functions/antifind.m
new file mode 100644
index 0000000..109b676
--- /dev/null
+++ b/functions/antifind.m
@@ -0,0 +1,29 @@
+%% Converts numerical indices to logical indices (opposite of built-in find)
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+function idxBool = antifind(idx, len)
+ idxBool = zeros(len, 1, 'logical');
+ idxBool(idx) = true;
+end
\ No newline at end of file
diff --git a/functions/crossChCorr.m b/functions/crossChCorr.m
new file mode 100644
index 0000000..fad3698
--- /dev/null
+++ b/functions/crossChCorr.m
@@ -0,0 +1,105 @@
+function rsqs = crossChCorr(tt, V, srate)
+%
+% This function calculates how mean cross-channel R-squared change at each CAR size
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+ % Ensure ch x time points structure if k = 1
+ if ismatrix(V)
+ if size(V, 1) == length(tt), V = V'; end
+ end
+
+ % Estimate sampling frequency from time points if necessary (not recommended)
+ if nargin < 3 || isempty(srate)
+ srate = (length(tt) - 1) / (max(tt) - min(tt));
+ warning('sampling frequency was estimated to be %0.02f', srate);
+ end
+
+ nChs = size(V, 1);
+ nTrs = size(V, 3);
+ assert(size(V, 2) == length(tt), 'Error: second dimension does not match tt. Data should be channels x timepoints x trials');
+
+ stats = struct();
+
+ % Create 60, 120, 180 Hz notch filters
+ for ff = 60:60:180
+ dNotch(ff/60) = designfilt('bandstopiir', 'FilterOrder', 4, ...
+ 'DesignMethod', 'butter', ...
+ 'HalfPowerFrequency1', ff-2, ...
+ 'HalfPowerFrequency2', ff+2, ...
+ 'SampleRate', srate);
+ end
+
+ % notch filter signals to use for covariance calculating
+ Vclean = V;
+ for ff = 60:60:180
+ for ii = 1:nTrs % filter per trial
+ Vclean(:, :, ii) = filtfilt(dNotch(ff/60), Vclean(:, :, ii)')';
+ end
+ end
+
+ % extract data on response interval to perform analysis on
+ twin = [0.01, 0.3]; % should this be an input parameter? For now keep hardcoded
+ Vseg = Vclean(:, tt >= twin(1) & tt <= twin(2), :);
+ VorigSeg = V(:, tt >= twin(1) & tt <= twin(2), :);
+
+ % Rank channels. If single trial rank based on increasing variance. If more than 1 trial, rank by covariance (self-self excluded). Better estimator of consistent structure
+ stats.vars = nan(nChs, 1);
+ for ii = 1:nChs
+ covCurr = cov(squeeze(Vseg(ii, :, :)));
+ stats.vars(ii) = mean(covCurr(logical(triu(ones(size(covCurr)), 1))), 'all');
+ end
+ [~, stats.order] = sort(stats.vars);
+
+ %
+ rsqs = nan(nChs, nChs, nChs+1); % first column corresponds to no rereferencing
+
+ % Pull subset U (order(1:ii)) from V, growing U each time by 1, to construct CAR
+ for ii = 0:nChs % min CAR size is 2
+
+ % rereferenced signal is calculated from the non-notch data
+ if ii == 0
+ VsegReref = VorigSeg;
+ else
+ VsegReref = VorigSeg - mean(VorigSeg(stats.order(1:ii), :, :), 1);
+ end
+
+ % apply notch filter to the rereferenced signal before calculating correlations, to reduce effect of line noise
+ for ff = 60:60:180
+ for jj = 1:nTrs % filter per trial
+ VsegReref(:, :, jj) = filtfilt(dNotch(ff/60), VsegReref(:, :, jj)')';
+ end
+ end
+
+ VsegRerefMean = mean(VsegReref, 3); % mean across trials
+
+ r = corr(VsegRerefMean');
+ r(1:(nChs+1):end) = nan; % omit self-self correlations
+
+ rsqs(:, :, ii+1) = r.^2; % convert to R^2 value
+
+
+ end
+
+end
diff --git a/functions/genRandSig.m b/functions/genRandSig.m
new file mode 100644
index 0000000..c5126c5
--- /dev/null
+++ b/functions/genRandSig.m
@@ -0,0 +1,52 @@
+%% File to simulate an evoked potential as part of the code to test CARLA
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+function sig = genRandSig(tt, nChs, amp)
+ % tt = time points
+ % amp = average amplitude of s1
+ % nChs = number of channels to generate
+ % Adds s1 and s2. s1 represents faster, direct response. s2 represents slower, indirect response
+
+ % signal parameters
+ A = 0.8*amp + 0.4*amp*rand(nChs, 1); % amplitudes
+
+ tau1 = 0.01 + rand(nChs, 1)*0.02; % time constant to dampen the faster sinusoid (s1)
+ tau2 = 0.06 + rand(nChs, 1)*0.08;
+ tauOn1 = 0.005; % time constant of subtracted exponential decay for both
+ tauOn2 = 0.025; % time constant of subtracted exponential decay for both
+
+ f1 = 8 + rand(nChs, 1)*4; % freq (hz) of first sinusoid (mean period = 100 ms)
+ f2 = 1 + rand(nChs, 1)*2; % freq of second sinusoid, (mean period = 500 ms)
+
+ ph1 = rand(nChs, 1)*2*pi; % phase of first sinusoid
+ ph2 = rand(nChs, 1)*2*pi; % phase of second sinusoid
+
+ sig = zeros(length(tt), nChs);
+ for ii = 1:nChs
+ sig(:, ii) = A(ii)*( (exp(-tt/tau1(ii)) - exp(-tt/tauOn1)).*sin(2*pi*f1(ii)*tt - ph1(ii)) + (exp(-tt/tau2(ii)) - exp(-tt/tauOn2)).*sin(2*pi*f2(ii)*tt - ph2(ii)) );
+ end
+ sig(tt < 0, :) = 0; % make causal
+
+end
\ No newline at end of file
diff --git a/functions/genRandSigOrig.m b/functions/genRandSigOrig.m
new file mode 100644
index 0000000..a8603a6
--- /dev/null
+++ b/functions/genRandSigOrig.m
@@ -0,0 +1,51 @@
+%% Older, alternative version of generating random signals. Used for schematic purposes, not for quantification
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+function sig = genRandSigOrig(tt, nChs, amp)
+ % tt = time points
+ % amp = average amplitude of s1
+ % nChs = number of channels to generate
+ % Adds s1 and s2. s1 represents faster, direct response. s2 represents slower, indirect response
+
+ % signal parameters
+ A = 0.8*amp + 0.4*amp*rand(nChs, 1); % amplitudes
+
+ tau1 = 0.05 + rand(nChs, 1)*0.06; % time constant of slower exp decay
+ tau2 = ones(nChs, 1)*0.01; % time constant of second exp decay (shorter, fixed)
+ tau3 = 0.025; % time constant to dampen the nested, faster sinusoid (s1), corresponds to half the half-period of
+
+ f1 = 8 + rand(nChs, 1)*4; % freq (hz) of first sinusoid (mean period = 100 ms)
+ f2 = 2 + rand(nChs, 1)*2; % freq of second sinusoid, (mean period = 250 ms)
+
+ ph1 = rand(nChs, 1)*2*pi; % phase of first sinusoid
+ ph2 = rand(nChs, 1)*2*pi; % phase of second sinusoid
+
+ sig = zeros(length(tt), nChs);
+ for ii = 1:nChs
+ sig(:, ii) = A(ii)*( (exp(-tt/tau1(ii)) - exp(-tt/tau2(ii))) .* (exp(-tt/tau3).*sin(2*pi*f1(ii)*tt - ph1(ii)) + sin(2*pi*f2(ii)*tt - ph2(ii))));
+ end
+ sig(tt < 0, :) = 0; % make causal
+
+end
\ No newline at end of file
diff --git a/functions/getNeighborChs.m b/functions/getNeighborChs.m
new file mode 100644
index 0000000..e576fae
--- /dev/null
+++ b/functions/getNeighborChs.m
@@ -0,0 +1,57 @@
+%% Returns the neighbors within n contacts for input channels
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+function neighbors = getNeighborChs(chs, n, selfincl)
+ if nargin < 3, selfincl = true; end % whether or not the input channel(s) should be included as neighbors
+ if nargin < 2 || isempty(n), n = 1; end
+
+ isDigit = @(x) x > 47 & x < 58; % returns true for char array elements that are digits (0 - 9)
+
+ if ischar(chs), chs = {chs}; end % if only 1 channel is requested
+
+ neighbors = {};
+ for ii = 1:length(chs)
+ ch = strip(upper(chs{ii}));
+
+ % split channel name into lead and contact
+ lead = ch(~isDigit(ch));
+ contact = str2double(ch(isDigit(ch)));
+
+ % all neighoring contact integers within n
+ neighborcontacts = contact + (-n:n);
+ neighborcontacts(neighborcontacts < 1 | neighborcontacts > 18) = [];
+
+ neighbors = [neighbors, arrayfun(@(x) sprintf('%s%d', lead, x), neighborcontacts, 'UniformOutput', false)];
+ end
+
+ % remove redundant entries
+ neighbors = unique(neighbors, 'stable');
+
+ % remove input channels if necessary
+ if ~selfincl
+ neighbors = setdiff(neighbors, chs);
+ end
+
+end
\ No newline at end of file
diff --git a/functions/rmBadTrialsAnnots.m b/functions/rmBadTrialsAnnots.m
new file mode 100644
index 0000000..ad50720
--- /dev/null
+++ b/functions/rmBadTrialsAnnots.m
@@ -0,0 +1,57 @@
+function dataOut = rmBadTrialsAnnots(dataIn, chs, annots)
+%% Removes bad trials according to text in a cell array of chars, annotations
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+ assert(size(dataIn, 1) == length(chs), 'Error: First dimension of input data does not match channel names');
+ assert(size(dataIn, 3) == length(annots), 'Error: Third dimension of input data does not match annotations length');
+
+ dataOut = dataIn;
+
+ annots = upper(strip(annots)); % clean up, make uppercase
+
+ % Remove trials that are bad across all channels
+ trs2Rm = strcmp(annots, 'ALL');
+ fprintf('Removing %d/%d trials\n', sum(trs2Rm), length(annots));
+ dataOut(:, :, trs2Rm) = [];
+ annots(trs2Rm) = [];
+
+ % find channels bad in individual trials and set as nan, to be removed later
+ for ii = 1:length(annots)
+ try
+ annot = annots{ii};
+
+ % Channels mentioned in description. matches strs starting with L/R, followed by 1 or more A-Z letters, then 1 or more numerical digits
+ out = regexp(annot, '[LR][A-Z]+\d+', 'match');
+ if isempty(out), continue; end
+
+ % identify which channels are bad at this trial
+ idxbad = ismember(chs, out);
+ fprintf('Setting %d bad channels at trial number %d to NaN values: %s\n', sum(idxbad), ii, string(join(out, ', ')));
+ dataOut(idxbad, :, ii) = nan;
+ catch
+ end
+ end
+
+end
\ No newline at end of file
diff --git a/main.m b/main.m
new file mode 100644
index 0000000..5253c9a
--- /dev/null
+++ b/main.m
@@ -0,0 +1,179 @@
+%% Use each section of this main script to output all results and figures for the manuscript.
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package: main script
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Unless specified, each section of code below can be run directly from this script to generate all outputs.
+% The scripts listed may also be opened so that individual sections within can be separately examined and run.
+% Note that in output filenames, "zminmean" refers to zeta in the manuscript
+%
+% Subject CCEP and anatomical data will be available in BIDS format upon manuscript publication.
+% In the meantime, all results on simulated data can be reproduced with the appropriate section below.
+%
+%% Configure paths and dependencies
+
+% Make sure you are currently working in the root manuscript directory (CARLA)
+addpath('functions');
+
+% Data should be downloaded from OpenNeuro and copied to a folder here (./data)
+dataPath = 'data';
+
+% Add paths to other dependencies. Change the dummy paths listed here to where they are located on your system.
+addpath(genpath('path/to/mnl_ieegBasics'));
+addpath(genpath('path/to/matmef'));
+addpath('path/to/mnl_seegview');
+addpath(genpath('path/to/vistasoft'));
+addpath('path/to/spm12'); addpath(fullfile('path/to/spm12', 'config'));
+addpath(genpath('path/to/freesurfer/7.1.1/matlab'));
+addpath(genpath('path/to/gifti'));
+
+% To save figures as vectorized objects
+set(0, 'DefaultFigureRenderer', 'painters');
+
+%% i) Generate T1 MRI slices with stimulation site electrodes overlaid (Figures 4A, 6C)
+
+% This was done using sliceGUI.m in the mnl_seegview repository:
+% https://github.com/MultimodalNeuroimagingLab/mnl_seegview
+
+% Define paths to defaced T1 MRI and matching electrodes table
+niiPath = fullfile(dataPath, 'derivatives', 'freesurfer', 'sub-1', 'sub-1_ses-mri01_T1w_acpc_deFaced.nii');
+electrodes = fullfile(dataPath, 'sub-1', 'ses-ieeg01', 'ieeg', 'sub-1_ses-ieeg01_electrodes.tsv');
+
+% Run sliceGUI, and manually interact with it
+sliceGUI(niiPath, electrodes);
+
+%% Table 1: Calculate the number of measurement electrodes (channels) used in each subject
+
+% outputs are directly printed to command window
+tab1_dispNumElectrodes
+
+%% Figure 1: Save simulated CCEP channels and their theoretical components
+
+% outputs saved to .\output\simComponents
+fig1_simCCEPComponents
+
+%% Figure 2: Methods illustration of how CARLA works, on single-trial simulated CCEP data
+
+% outputs saved to .\output\simulationTrial
+fig2_CARLAonSingleTrial
+
+%% Figure 3: Save example of enveloped sinusoids used to simulate evoked potentials
+
+% outputs saved to .\output\EPConstruct
+fig3_EPConstruct
+
+%% Figure 4: Simulate 50 channels with responsiveness varying from 0 to 50, with 30 repetitions each, and apply CARLA
+
+% First run this to generate outputs for each number of responsive channels
+% outputs saved to .\output\simLoop
+Aglobal = 0; % no global signal desired (0 amplitude)
+simCCEPLoop;
+
+% Then, keeping Aglobal in the workspace, run this to make summary figures (4C, D), as well as identify median candidates for 4A, B
+% .\output\simLoop will be renamed .\output\simLoopNoGlob, where outputs are saved.
+summarizeSimCCEPLoop
+
+%% Figure 5: Apply CARLA to real data (sub 1, stim site 1)
+
+% Configure the subject, data run, and stim site
+sub = 'sub-1';
+run = 2;
+site = 'RMO8-RMO9';
+
+% which channels to plot in panel F
+chs2Plot = [100, 203];
+
+% outputs saved to .\output\realCCEPs\. These correspond to figures 5B-F. To generate 5A (stim site overlaid on T1 MRI), see section i), above.
+applyCARLARealCCEPs
+
+%% Figure S1: Average channel responses for sub 1, stim site 1 after applying CAR vs. CARLA
+
+% This supplemental figure is a direct addendum to Figure 5, showing the mean signal for the first 50 channels after rereferencing by CARLA or CAR
+% Run after the previous section, keeping all workspace variables
+figS1_compareCARvCARLA
+
+%% Figure S2: Apply CARLA to real data at another 4 stim sites in subjects 1-4
+
+% config parameters for each subject shown in S1
+config = struct('sub', {'sub-1', 'sub-2', 'sub-3', 'sub-4'}, ...
+ 'run', {1, 12, 3, 2}, ...
+ 'site', {'ROC11-ROC12', 'RZ4-RZ5', 'LA2-LA3', 'RO3-RO4'});
+chs2Plot = 1; % dummy value, no individual channels are plotted in Figure S1
+
+for ii = 1:length(config)
+ sub = config(ii).sub;
+ run = config(ii).run;
+ site = config(ii).site;
+
+ % outputs saved to .\output\realCCEPs\.
+ % Omission of the middle section of sorted channels (subjects 1-3) in Figure S1 was done manually in illustrator.
+ applyCARLARealCCEPs
+end
+
+%% Figure S3: Save simulated CCEP channels, this time with a global signal, and their theoretical components
+
+% outputs saved to .\output\simComponentsGlobal
+figS3_simCCEPComponentsGlobalSig
+
+%% Figure 6A, B Simulate 50 channels with responsiveness from 0 to 50, like Figure 4, but now with a global signal added
+
+% First run this to generate outputs for each number of responsive channels
+% outputs saved to .\output\simLoop
+Aglobal = 25; % global signal with mean amplitude 25
+simCCEPLoop;
+
+% Then, keeping Aglobal in the workspace, run this to make summary figure (6B), and to identify a trial candidate for 6A
+% .\output\simLoop will be renamed .\output\simLoopWithGlob, where outputs are saved.
+summarizeSimCCEPLoop
+
+% To generate figure 6C, see section i) on using sliceGUI, above.
+
+%% Figure 6D-F: Apply CARLA to real data containing global signal (sub 1, stim site 2)
+
+% Configure the subject, data run, and stim site
+sub = 'sub-1';
+run = 3;
+site = 'RK5-RK6';
+
+% which channels to plot in panel F
+chs2Plot = [23, 51, 205];
+
+% outputs saved to .\output\realCCEPs\
+applyCARLARealCCEPs
+
+%% Figure 7: Loop through stimulation sites in each subject and compute optimal CARLA cutoff
+% This also calculates the cross-channel R-squared values for all possible adjusted common average sizes, to be summarized in Figure 8
+
+% First run this script. Outputs are saved to .\output\realCCEPsLoop\.
+applyCARLARealCCEPsLoop
+
+% Then run this script to generate outputs for Figure 7
+fig7_summarizeCARLARealCCEPs
+
+%% Figure 8: Calculate cross-channel R-squared values for no CAR vs different types of fixed CARs vs CARLA
+
+% outputs saved to .\output\realCCEPsLoop\.
+fig8_summarizeInterChCorr
+
diff --git a/simCCEPLoop.m b/simCCEPLoop.m
new file mode 100644
index 0000000..adffd52
--- /dev/null
+++ b/simCCEPLoop.m
@@ -0,0 +1,175 @@
+%% This script creates N channels with a certain number of responsive channels.
+% At each loop it tests results and accuracy of CARLA.
+% Set Aglobal to 0 for outputs to Figure 4 and set Aglobal to 25 for outputs to Figure 6A,B
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+%% First control the parameters of the simulated channels
+
+srate = 4800;
+tt = -0.5:1/srate:1-1/srate;
+
+nchs = 50; % number of simulated channels
+ntrs = 12; % number of trials
+
+cmSens = [1, 165/255, 0]; % use orange as color for the sensitive optimum
+
+%% Loop through number of responsive channels
+
+rng('default');
+
+for nresp = 0:1:50
+
+ %% Create output dir for current responsiveness
+
+ fprintf('%d of %d channels responsive\n', nresp, nchs);
+
+ outdir = fullfile('output', 'simLoop', sprintf('nchs%d-%d', nchs, nresp));
+ mkdir(outdir);
+
+ chsResp = 1:nresp; % first x channels responsive, for ease. order doesn't matter. will sort and color them red
+
+ reps = 30;
+
+ %%
+ for rr = 1:reps % multiple repetitions at each responsive level
+
+ fprintf('.'); % loading bar
+
+ % i) artifact only
+ V0 = zeros(length(tt), nchs);
+ Aart = 50 + rand(nchs, 1)*5; % slightly different artifact amplitudes for each channel
+ artifact = sin(2*pi*600*tt)';
+ artifact(tt < 0 | tt > 0.002) = 0;
+ V0 = V0 + artifact*Aart';
+
+ % A) Add the evoked potentials
+ A = 100;
+ V1 = V0;
+ sig = genRandSig(tt, length(chsResp), A);
+ V1(:, chsResp) = V0(:, chsResp) + sig;
+
+ % B) Option to add a global noise.
+ % Configure in main script for Aglobal to be 0 (for figure 4) or 25 (for figure 6A-B)
+ if Aglobal
+ sigCommon = genRandSig(tt, 1, Aglobal);
+ else
+ sigCommon = 0; % no global signal
+ end
+ V2 = V1 + sigCommon;
+
+ % C) Add common noise to all channels at each trial
+ V3 = repmat(V2', 1, 1, ntrs); % ch x time points x trial
+ phLN = rand(ntrs, 3)*2*pi; % LN phases
+ LN = zeros(length(tt), ntrs);
+ for ii = 1:ntrs
+ LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3));
+ end
+
+ % brown noise component of common noise shared across channels
+ BN = cumsum(0.4*randn(2*length(tt), ntrs));
+ BN = ieeg_highpass(BN, srate, true);
+ BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :);
+
+ noiseCommon = LN + BN;
+ V3 = V3 + shiftdim(noiseCommon, -1);
+
+ % D) add random brown noise
+ noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it
+ for ii = 1:nchs
+ noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true);
+ end
+ noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :);
+ V4 = V3 + noiseRand;
+
+ %% Apply CARLA and plot outputs
+
+ [Vout, CAR, stats] = CARLA(tt, V4, srate, true); % get the sensitive output
+
+ % number of channels used for the CAR
+ nCAR = length(stats.chsUsed);
+ [~, nCARglob] = max(mean(stats.zMinMean, 2)); % number of channels at global maximum
+
+ % Plot average zmin across trials
+ figure('Position', [200, 200, 400, 300]); hold on
+ errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k-o');
+ plot(nCARglob, mean(stats.zMinMean(nCARglob, :), 2), 'b*'); % global max as blue
+ if nCARglob ~= nCAR; plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens); end
+ yline(0, 'Color', 'k');
+ saveas(gcf, fullfile(outdir, sprintf('zmin_%d-%d_rep%d', nchs, nresp, rr)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('zmin_%d-%d_rep%d', nchs, nresp, rr)), 'svg');
+
+ % Plot the variance of channels, sorted in increasing order
+ vars = stats.vars(stats.order);
+ figure('Position', [200, 200, 400, 300]); hold on
+ plot(vars, 'k-o', 'LineWidth', 1, 'MarkerFaceColor', 'k');
+ xline(nCARglob + 0.5, 'Color', 'b'); % global threshold
+ if nCARglob ~= nCAR; xline(nCAR + 0.5, 'Color', cmSens); end % sensitive threshold
+ xlim([0, nchs+1]);
+ saveas(gcf, fullfile(outdir, sprintf('covar_%d-%d_rep%d', nchs, nresp, rr)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('covar_%d-%d_rep%d', nchs, nresp, rr)), 'svg');
+
+ % Sort and plot channels by increasing covariance, draw line at cutoff
+ V4MeanSorted = mean(V4(stats.order, :, :), 3);
+
+ % create colormap where responsive channels are red
+ respBool = antifind(chsResp, nchs);
+ respBool = respBool(stats.order); % logical array of where responsive channels are
+ cm = zeros(nchs, 3);
+ cm(respBool, 1) = 1; % make red
+ figure('Position', [200, 200, 250, 600]);
+ yspace = 80;
+ ys = ieeg_plotTrials(tt, V4MeanSorted', 80, [], cm, 'LineWidth', 1);
+ yline(ys(nCARglob)-yspace/2, 'Color', 'b', 'LineWidth', 1.5);
+ if nCARglob ~= nCAR; yline(ys(nCAR)-yspace/2, 'Color', cmSens, 'LineWidth', 1.5); end
+ xlim([-0.1, 0.5]); set(gca, 'xtick', [0, 0.5]);
+ xlabel('Time (s)'); ylabel('Channels');
+ saveas(gcf, fullfile(outdir, sprintf('chsSorted_%d-%d_rep%d', nchs, nresp, rr)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('chsSorted_%d-%d_rep%d', nchs, nresp, rr)), 'svg');
+
+ close all;
+
+ % Accuracy values. positive means responsive/excluded from CAR
+ TP = sum(find(respBool) > nCAR); % responsive channels successfully excluded from CAR (above the cutoff)
+ TN = sum(find(~respBool) <= nCAR); % NR channels successfully below or at cutoff
+ FN = sum(find(respBool) <= nCAR); % responsive channels incorrectly included in CAR. *This matters most
+ FP = sum(find(~respBool) > nCAR); % NR channels incorrectly excluded from CAR
+
+ % same for the global threshold
+ TPglob = sum(find(respBool) > nCARglob);
+ TNglob = sum(find(~respBool) <= nCARglob);
+ FNglob = sum(find(respBool) <= nCARglob);
+ FPglob = sum(find(~respBool) > nCARglob);
+
+ fid = fopen(fullfile(outdir, sprintf('accuracy_%d-%d_rep%d.txt', nchs, nresp, rr)), 'w');
+ fprintf(fid, 'TP\t%d\nTN\t%d\nFN\t%d\nFP\t%d\n', TP, TN, FN, FP);
+ fprintf(fid, 'TPglob\t%d\nTNglob\t%d\nFNglob\t%d\nFPglob\t%d', TPglob, TNglob, FNglob, FPglob);
+ fclose(fid);
+
+ end
+
+ fprintf('\n');
+
+end
+
diff --git a/summarizeSimCCEPLoop.m b/summarizeSimCCEPLoop.m
new file mode 100644
index 0000000..c3ee41c
--- /dev/null
+++ b/summarizeSimCCEPLoop.m
@@ -0,0 +1,197 @@
+%% This script takes outputs from simCCEPLoop, and evaluates accuracy at each level of responsiveness.
+% Generates outputs for Figure 4C,D and 6B
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+
+%% Configure loading directories
+%
+nchs = 50; % total number of channels
+
+% rename the output folder depending on whether Aglobal was 0 to avoid overwriting Figure 4 outputs with Figure 6 outputs
+if ~Aglobal % no global signal
+ rootdir = fullfile('output', 'simLoopNoGlob');
+else % yes global signal
+ rootdir = fullfile('output', 'simLoopWithGlob');
+end
+movefile(fullfile('output', 'simLoop'), rootdir);
+
+%% Read and store accuracy values at each subset of responsiveness
+
+% each field is reps x nresp (starting at 0). XXglob = at cutoff found using global maximum
+accuracy = struct();
+accuracy.TP = nan(30, nchs+1);
+accuracy.TN = nan(30, nchs+1);
+accuracy.FN = nan(30, nchs+1);
+accuracy.FP = nan(30, nchs+1);
+accuracy.TPglob = nan(30, nchs+1);
+accuracy.TNglob = nan(30, nchs+1);
+accuracy.FNglob = nan(30, nchs+1);
+accuracy.FPglob = nan(30, nchs+1);
+
+for nresp = 0:nchs
+
+ fprintf('.');
+
+ indir = fullfile(rootdir, sprintf('nchs%d-%d', nchs, nresp));
+
+ for rep = 1:30
+
+ T_acc = readtable(fullfile(indir, sprintf('accuracy_%d-%d_rep%d.txt', nchs, nresp, rep)), 'FileType', 'text', 'Delimiter', '\t');
+ accuracy.TP(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TP'));
+ accuracy.TN(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TN'));
+ accuracy.FN(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FN'));
+ accuracy.FP(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FP'));
+ accuracy.TPglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TPglob'));
+ accuracy.TNglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TNglob'));
+ accuracy.FNglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FNglob'));
+ accuracy.FPglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FPglob'));
+
+ end
+
+end
+fprintf('\n');
+
+%% Identify examples that correspond to median values at each # responsive channels
+
+FNmed = median(accuracy.FN);
+FPmed = median(accuracy.FP);
+FNglobMed = median(accuracy.FNglob);
+FPglobMed = median(accuracy.FPglob);
+
+if contains(rootdir, 'simLoopWithGlob')
+ idxes = 26; % 50% for main figure with global signal
+ FNmedSelect = FNmed(idxes);
+ FPmedSelect = FPmed(idxes);
+ FNglobmedSelect = FNglobMed(idxes);
+ FPglobmedSelect = FPglobMed(idxes);
+ fprintf('Median FN at 25 resp chs = %0.1f\n', FNmedSelect);
+ fprintf('Median FP at 25 = %0.1f\n', FPmedSelect);
+ fprintf('Median FN (global) at 25 = %0.1f\n', FNglobmedSelect);
+ fprintf('Median FP (global) at 25 = %0.1f\n', FPglobmedSelect);
+
+elseif contains(rootdir, 'simLoopNoGlob')
+ idxes = [1, 11, 21, 31, 41, 46]; % every 20%, for main figure with no global signal
+ FNmedSelect = FNmed(idxes);
+ FPmedSelect = FPmed(idxes);
+ FNglobmedSelect = FNglobMed(idxes);
+ FPglobmedSelect = FPglobMed(idxes);
+ fprintf('Median FN at 0, 10, 20, 30, 40, 45 resp chs = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FNmedSelect(:));
+ fprintf('Median FP at 0, 10, 20, 30, 40, 45 = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FPmedSelect(:));
+ fprintf('Median FN (global) at 0, 10, 20, 30, 40, 45 = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FNglobmedSelect(:));
+ fprintf('Median FP (global) at 0, 10, 20, 30, 40, 45 = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FPglobmedSelect(:));
+
+end
+
+% which rep to use that best matches the median
+repSelect = zeros(size(idxes));
+for ii = 1:length(idxes)
+ rep = find(accuracy.FN(:, idxes(ii)) == floor(FNmedSelect(ii)) & ...
+ accuracy.FP(:, idxes(ii)) == floor(FPmedSelect(ii)) & ...
+ accuracy.FNglob(:, idxes(ii)) == floor(FNglobmedSelect(ii)) & ...
+ accuracy.FPglob(:, idxes(ii)) == floor(FPglobmedSelect(ii)), 1, 'first');
+ if isempty(rep), rep = nan; end
+ repSelect(ii) = rep;
+end
+fprintf('Rep to use at responsiveness: %d \n', repSelect(:));
+
+%% Absolute values of FP, FN
+
+cmSens = [1, 165/255, 0]; % orange = sensitive threshold
+
+% Using global optimum
+% up to 90% responsiveness (45/50 responsive channels)
+figure('Position', [200, 200, 600, 600]);
+subplot(2, 1, 1); % FN = responsive channels included in the CAR
+boxplot(accuracy.FNglob(:, 1:46), 'Positions', 0:45, 'Colors', 'b', 'PlotStyle', 'compact');
+xlim([-1, 46]); ylim([-1, 50]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off');
+ylabel('Number of incorrectly assigned channels');
+title('False Negative (Resp. Chs in CAR)');
+
+subplot(2, 1, 2);
+boxplot(accuracy.FPglob(:, 1:46), 'Positions', 0:45, 'Colors', 'b', 'PlotStyle', 'compact');
+xlim([-1, 46]); ylim([-1, 50]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off');
+xlabel('Number of responsive channels / 50');
+ylabel('Number of incorrectly assigned channels');
+title('False Positive (Non-Resp. Chs in CAR)');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFPglob', nchs)), 'png');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFPglob', nchs)), 'svg');
+
+% For sensitive method
+% up to 90% responsiveness (45/50 responsive channels)
+figure('Position', [200, 200, 600, 600]);
+subplot(2, 1, 1); % FN = responsive channels included in the CAR
+boxplot(accuracy.FN(:, 1:46), 'Positions', 0:45, 'Colors', cmSens, 'PlotStyle', 'compact');
+xlim([-1, 46]); ylim([-1, 50]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off');
+ylabel('Number of incorrectly assigned channels');
+title('False Negative (Resp. Chs in CAR)');
+
+subplot(2, 1, 2);
+boxplot(accuracy.FP(:, 1:46), 'Positions', 0:45, 'Colors', cmSens, 'PlotStyle', 'compact');
+xlim([-1, 46]); ylim([-1, 50]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off');
+xlabel('Number of responsive channels / 50');
+ylabel('Number of incorrectly assigned channels');
+title('False Positive (Non-Resp. Chs in CAR)');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFP', nchs)), 'png');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFP', nchs)), 'svg');
+
+%% Sensitivity specificity plot
+
+% sensitivity and specificity
+sensglob = accuracy.TPglob ./ (accuracy.TPglob + accuracy.FNglob);
+specglob = accuracy.TNglob ./ (accuracy.TNglob + accuracy.FPglob);
+sens = accuracy.TP ./ (accuracy.TP + accuracy.FN); % more important. need to sensitively detect responsive channels
+spec = accuracy.TN ./ (accuracy.TN + accuracy.FP);
+
+% balanced accuracy
+accglob = (sensglob + specglob)/2;
+acc = (sens + spec)/2;
+
+
+% Using global optimum
+cmBlue = [0, 0, 139; 115, 147, 179]/255;
+figure('Position', [200, 200, 600, 300]); hold on
+errorbar(0:45, mean(sensglob(:, 1:46)), std(sensglob(:, 1:46)), '-o', 'MarkerSize', 5, 'MarkerFaceColor', cmBlue(1, :), 'Color', cmBlue(1, :));
+errorbar(0:45, mean(specglob(:, 1:46)), std(specglob(:, 1:46)), '--s', 'MarkerSize', 6, 'Color', cmBlue(2, :));
+xlim([-1, 46]); ylim([-0.2, 1.2]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50);
+xlabel('Number of responsive channels / 50');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpecglob', nchs)), 'png');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpecglob', nchs)), 'svg');
+
+% Using sensitive method
+cmOrange = [139, 64, 0; 255, 172, 28]/255;
+figure('Position', [200, 200, 600, 300]); hold on
+errorbar(0:45, mean(sens(:, 1:46)), std(sens(:, 1:46)), '-o', 'MarkerSize', 5, 'MarkerFaceColor', cmOrange(1, :), 'Color', cmOrange(1, :));
+errorbar(0:45, mean(spec(:, 1:46)), std(spec(:, 1:46)), '--s', 'MarkerSize', 6, 'Color', cmOrange(2, :));
+xlim([-1, 46]); ylim([-0.2, 1.2]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50);
+xlabel('Number of responsive channels / 50');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpec', nchs)), 'png');
+saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpec', nchs)), 'svg');
+
diff --git a/tab1_dispNumElectrodes.m b/tab1_dispNumElectrodes.m
new file mode 100644
index 0000000..f04c50e
--- /dev/null
+++ b/tab1_dispNumElectrodes.m
@@ -0,0 +1,55 @@
+%% This script calculates and prints the number of channels used in the analysis for each subject ("good" in at least 1 run of data)
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 Harvey Huang
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see .
+%
+% path to subject data
+dataPath = 'data';
+
+% subjects and runs corresponding to each subject to load
+subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'};
+runs = {[1, 2, 3, 5], [12, 13, 14], [2, 3, 4], [1, 2]};
+
+% iterate through each subject, then each run, and calculate number of good channels
+for ss = 1:4
+
+ sub = subs{ss};
+
+ seegChStatus = {};
+
+ % iterate across runs for each subject
+ for rr = 1:length(runs{ss})
+
+ run = runs{ss}(rr);
+
+ channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run));
+ channels = ieeg_readtableRmHyphens(channelsPath);
+ seegChStatus = [seegChStatus, channels.status(strcmp(channels.type, 'SEEG'))];
+
+ end
+
+ % Good channels are defined as any that are not bad across all runs (and therefore not used in any analysis)
+ goodChs = any(strcmp(seegChStatus, 'good'), 2);
+ fprintf('Number of good sEEG channels across all runs in %s = %d\n', sub, sum(goodChs));
+
+end