preproc_v4

# MAIN for expt2 - Dec 30, 2012
# This main file will eventually be the core of data analysis pipeline
echo
echo
echo SWDfinder Data Pre-processing
echo ~f~ version 1.0 created December 31, 2012
echo 
# Copyright (C) 2014 Christian Richard

# 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 <http://www.gnu.org/licenses/gpl.html>. 

# Author contact info:
# lifepupil@gmail.com

# - must cross check values in {ref tannot files} with their expt1 counterparts to make sure they are identical,
# - investigate Channel Index in annotation file to see whether they correspond to channel where SWD was noted by Wayne, Verity, &c.
#	if yes, then these values can be used to increase specificity of the 'best' SWDs from which to do reverse engineering type project.
# - make sure that all required function calls are made within each pre-processing step
# - change code that uses the .pd1m files (m=2,3,4) AND resave expt1 psync folder files with corrected m's
# - make GUI


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# INSTRUCTIONS FROM HARMONIE TO LASTWAVE
# This function to be used on Harmonie data output generated by:
#	1. In Harmonie, choose "patient"
#	2. Right mouse click on desired file and select Properties
# 	3. Copy the mouse IDs and date (mmddyy) listed in the Identification field
# 	4. Make sure checkmark next to Workspaces | All Mice
# 	5. Select Tools | Section Marker
# 	6. Hold down left mouse button anywhere over montage and drag the line that appears to the leftmost point before releasing mouse button (total duration should be 0)
# 	7. Use slider at bottom of window to scroll to the end of the recording, then click anywhere on empty white part of screen
# 	8. Select File | Save Selection As Text... and paste contents of Identification field (step#3) into file name field
# 	9. Finally, to get annotation file, select File | Create Report then choose Event Distribution Workbook Report, then Clear SWDs and BIPOLAR. 
#	   Save the entire workbook using info copied in #3. Copy the worksheet with actual annotation values into a new book, then save as a .CSV file using same #3 info.
#	(OLD STEP - IGNORE) 10. Run HarmonieEEG_LW function to start
#	(OLD STEP - IGNORE) 11. Run annotConv to complete the annotation conversion
#	10. Move the EEG and annotation files (generated above) into the rawdata folder (for location check/set HarmonieReadPath in PATHS below).
#	11. Check main1 in a text editor to make sure that the flow controls are set to values reflecting the desired SWDfinder files.
#	12. Open lastwave and type "source main1" (without quotes) at the command prompt.
#	13. When finished, an unambiguous DONE message will appear in the lastwave window. Files will be found in folders as the PATHS set below.
#
# UPDATE (4-8-14):
# - broke commonPath into commonPath and commonFolder so that a more abbreviated identifier of strain can be passed to pd3d_v2
# - added PhaseDistSeq option to write phase difference coordinates of SWCs into text file for clustering algorithm in R
# - incorporated changes required for a fully realized PhaseDistSeq into pd3d_v2, including the creation of a trajectory velocity function trajVel
#
# UPDATE (5-6-14):
# - Modified PhaseDistSeq (and downstream functions pd3d_v2 and pd3d_v3.R) to write pertinent values at each timepoint to text file, including power at fundamental frequency
#	as requested by Alain and Benji
#
# TO DO (5-16-14):
# 1. ADD FEATURE TO COMPARE OVERLAP BETWEEN TWO PHASE DIFF DISTRIBUTIONS (SEE NOTES, row by row; A*B, take .X and count for coord %, add lowest val (A or B) to trajectory count)
# 
# UPDATE (5-30-14):
# - Enough modifications to justify saving as preproc_v2. 
# NOTABLY - changes to pd3d_v2 (now source pd3d_v3) permits generation of distribution whose coordinate values represent the NUMBER of SWCs that pass through it.
#	If a SWC had more than one timepoint with the same pd coordinates then the previous distribution method (-add) would would plot them all. The new
#	method (-add_by_swc) was developed so that distributions would more closely reflect the SWC content of, e.g. different strains (genetics) 
#	and EEG channels (brain regions).
# 
# UPDATE (7-25-14):
# - created pd3d_v4 to work with new preproc version (preproc_v4) on the PhaseDistSeq optSWDion
# - now a single text file containing SWD information at each timepoint (one per row, in order) for an entire strain
# - did all of this to permit more selective 3-D plotting with associated R script
# - creating a pd3d_v4.R to work with the above.
# UPDATE (7-31-14):
# - ADDED option to generate pd distribution of parts of SWDs, e.g. first 2 seconds of SWDs, to test whether the generally spikier SWCs at the start of a SWD 
# 	will have trajectories close to pure spike coordinates. This is done by changing manual annotation signal prior to calling pd3d function, and is activated
# by setting partialSWD to 1, and choosing sample range in partialSWDrng.
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# TIME KEEPING
startt = [time current -l]
if (startt[1]<10) {startt[1] = '0$startt[1]'} else {startt[1] = '$startt[1]'}
if (startt[2]<10) {startt[2] = '0$startt[2]'} else {startt[2] = '$startt[2]'}

#source postproc_v1

# >>~~~~~~~~~~~~~~~~~~~~>
# INSTANCE VARIABLES
os = 'win'
hours = 2
# commonFolder = {"expt2_filtered_LP03_HP50" "expt3_filtered_LP03_HP50" "Gabrg2_filtered_LP03_HP50" "Scn8a_filtered_LP03_HP50" "Scn8a_hom_filtered_LP03_HP50" "SMJ_filtered_LP03_HP50"}
commonFolder = {"Scn8a_filtered_LP03_HP50"}
freqIndexFile = 'freqind_amin2_vox10_omax6.txt'

# flow controls
createLWRs = 0
createAnnot = 0
scalGen = 0
refMaker = 0
phaseAnalysis = 0
PhaseDistHist = 0
	partialSWD = 1
	partialSWDrng = 0:200
	pd3dSWCs = 0
PhaseDistSeq = 0

# PhaseDistHist - generates phase difference distribution (pd3dSWCs setting requires that this process is activated)
# pd3dSWCs - controls whether all coordinates are plotted in -add method of pd3d_v3 (pd3dSWCs=0) or only each instance once by SWCs using -add_by_swc (pd3dSWCs=1)
# PhaseDistSeq - generates text files of SWC phase differences, power, etc. for use in R
# phaseCoordSeek - quick way to open existing distribution for display; also returns max values across distribution axis (currently phi1:2 only)

# ISI variables (from fundamental frequency phase)
maxVar = 0.15
minISInum = 2
maxISIdur = 40
minISIdur = 20

# Wavelet transform variables
sampleRate = 200
MorletScale = -1
amin = 2
vox = 10
omax = 6
morletOscillations = 3
MorletWindowSize = sampleRate/2

# harmonics analysis variables
# FOR SMJ 5.440941021 to 14.35872944 Hz
# startXcorr = 17
# endXcorr = 31
# FOR scn8a homozygotes 3.847326292 to 7.694652583 Hz
# startXcorr = 12
# endXcorr = 22
# ORIGINAL
startXcorr = 13
endXcorr = 26
# WIDEST RANGE
# startXcorr = 13
# endXcorr = 31

# MUST FINISH INTEGRATING THIS INTO CODE SO THAT DIFFERENT REF FILES CAN BE ACCESSED BY STRAIN
# ALSO INCLUDE AN IF STATEMENT TO TOGGLE BETWEEN OLD REF FILE NAMING SCHEME AND NEW ONE USING refParamStr STRING
# start_endXcorr = {13:26 13:26 13:26 13:26 12:22 17:31}

# 3-D phase difference histogram generation variables
# To selectively populate histogram by channel
frontbackRanges = <0,5>
diagRanges = <1,4>
#channelRanges = <frontbackRanges,diagRanges>
channelRanges = <0:5>
# /INSTANCE VARIABLES


# CONSTANTS
startX = 0
ISIdurations = 1
ranges = 0
durations = 1
SWDgroup = 0

sessionDuration = sampleRate * (hours*60*60)

# file name identifiers
WTparamStr = 'amin_'+'$amin'+'_vox_'+'$vox'+'_omax_'+'$omax'+'_mO_'+'$morletOscillations'+'_mscale_'+'$MorletScale'+'_winsize_'+'$MorletWindowSize'+'_sampling_'+'$sampleRate'+'_'
RoSparamStr = 'maxVar_'+'$maxVar'+'_minISInum_'+'$minISInum'+'_maxISIdur_'+'$maxISIdur'+'_minISIdur_'+'$minISIdur'+'_startXcorr_'+'$startXcorr'+'_endXcorr_'+'$endXcorr'+'_sampling_'+'$sampleRate'+'_'
refParamStr = '_startXcorr_'+'$startXcorr'+'_endXcorr_'+'$endXcorr'+'_'

# List of parameters and their indices used in Morlet wavelet transform of signal
WTparams = {}
ami = 0
vxi = 1
oxi = 2
mosci = 3
mscli = 4
sri = 5
scapi = 6
scppi = 7
wtfni = 8

# default values:
# WTparams[ami]		amin = 2
# WTparams[vxi]		vox = 10
# WTparams[oxi]		omax = 6
# WTparams[mosci]	morletOscillations = 3
# WTparams[mscli]	MorletScale = -1
# WTparams[sri]		sample rate = 200 Hz
# WTparams[scapi]	amplitude scalogram path
# WTparams[scppi] 	phase scalogram path
# WTparams[wtfni]	wavelet transform file name identifier

WTparams+= amin
WTparams+= vox
WTparams+= omax
WTparams+= morletOscillations
WTparams+= MorletScale
WTparams+= sampleRate
# WTparams+= scalAPath
# WTparams+= scalPPath
WTparams+= WTparamStr

phi_12 = 0
phi_13 = 1
phi_14 = 2

# for ref files - must match where values are stored in harmonicRefs2 function
xcov_max = 0
xcov_FF_indices = 1
xcorr_max = 2
xcorr_FF_indices = 3
# /CONSTANTS

# THIS PERMITS ANALYSIS OF MULTIPLE STRAINS
# foreach strain 0:0 {
foreach strain 0:commonFolder.size-1 {
	
	thisFolder = commonFolder[strain]

	# >>~~~~~~~~~~~~~~~~~~~~>
	# PATHS
	if (os=='win') {
		# commonPath = "C:/Users/lifepupil/LastWave/1Chris/Frankel/" + "$thisFolder" + "/"
		drivePath = "D:/"
		commonPath = drivePath + "$thisFolder" + "/"
	} elseif (os=='linux') {
		drivePath = "/media/Elements/"
		commonPath = drivePath + "$thisFolder" + "/"
	#	commonPath = "/home/chris/Desktop/1Chris/" + "$thisFolder" + "/"
	}
	HarmonieReadPath = commonPath+"rawdata/"
	annotPath = commonPath+"convAnnot/"
	lwrPath = commonPath+"lwr/"
	scalAPath = commonPath+"scalA/"
	scalPPath = commonPath+"scalP/"
	refPath = commonPath+"ref/"
	hampPath = commonPath+"hamp/"
	hphzPath = commonPath+"hphz/"
	psyncPath = commonPath+"psync/"
	pd3dPath = commonPath+"phaseWeight3D/"
	# >>~~~~~~~~~~~~~~~~~~~~>
	# /PATHS


	echo WORKING ON $thisFolder
	echo
	# Call to function that puts Harmonie output text files into image data structures so LW can work on them
	# Make sure that Harmonie output files stored in location specified by readPath

	if (createLWRs==1) {
		# FUNCTION CALLS
		source HarmonieEEG_LW3
		# /FUNCTION CALLS
		echo 
		echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		echo Extracting EEG data from Harmonie output files into LastWave image files
			
		HarmonieEEG_LW3 HarmonieReadPath lwrPath mouseIDsinFileName sessionDuration
		echo
	}

	if (createAnnot==1) {
		# FUNCTION CALLS
		source HarmonieAnnot_LW3
		# /FUNCTION CALLS
		echo 
		echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		echo Processing annotation table from Harmonie for use in LastWave
		HarmonieAnnot_LW HarmonieReadPath annotPath sessionDuration
		echo
	}

	if (scalGen==1) {
		# FUNCTION CALLS
		source scalogramGen2
		# /FUNCTION CALLS
		echo 
		echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		echo Calculating scalograms of EEG data
		fNames = [file list lwrPath+'*.lwr']
		foreach session 0:fNames.size-1 {
			mouseID = fNames[session]
			mouseID[[str substr mouseID ".lwr"]] := ''
			
			thisLWR = <;>
			iread thisLWR lwrPath+fNames[session]
			scalogramImage = Zero((WTparams[oxi]*WTparams[vxi]),thisLWR.ncol)
			scalogramGen a thisLWR scalogramImage WTparams mouseID

			clear a
			delete scalogramImage
			delete thisLWR
		}
	}

	if (refMaker==1) {
		echo 
		echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		echo Starting reference file maker
		# FUNCTION CALLS
		source harmonicRefs2
		#source numdur
		# /FUNCTION CALLS

		# determine EEG files' data to analyze
		# puts names of all LastWave EEG recording (.lwr) files into fNames
		fNames = [file list lwrPath+"*.lwr"]
		fNames = [listv sort fNames]
		fNum = fNames.size-1

		foreach session 0:fNum {
			rawEEG = <0;0>
			iread rawEEG lwrPath+fNames[session]
			sessionLength = rawEEG.ncol
			delete rawEEG
			
			ref = Zero(4,sessionLength)
			
			mouseFile = fNames[session]
			mouseFile[[str substr mouseFile '.lwr']] := ''

			foreach channel channelRanges {
				scalA = <0;0>
				echo
				echo Reading amplitude scalogram for $mouseFile channel $[channel+1]
				chfname = mouseFile+'_$WTparamStr'+'$channel'+'.scalA'
				iread scalA scalAPath+chfname
				
				echo Identifying best harmonics candidate for each timepoint
				harmonicRefs2 scalA ref startXcorr endXcorr
				delete scalA
				
				iwrite ref refPath+"$mouseFile"+"$refParamStr"+"$channel"+".ref"
				ref[:;:] :=0
			}
		}
	}

	if (phaseAnalysis==1) {
		echo 
		echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		echo Starting phase analysis
		# FUNCTION CALLS
		source harmonicVals2
		source pd_s
		# /FUNCTION CALLS

		fNames = [file list lwrPath+"*.lwr"]
		fNames = [listv sort fNames]
		fNum = fNames.size-1
		
		foreach session 0:fNum {
			
			rawEEG = <0;0>
			iread rawEEG lwrPath+fNames[session]
			sessionLength = rawEEG.ncol
			channelNum = rawEEG.nrow
			delete rawEEG
			
			mouseFile = fNames[session]
			mouseFile[[str substr mouseFile '.lwr']] := ''
			
			hamp = Zero(5,sessionLength)
			hphz = Zero(4,sessionLength)
			
			foreach channel 0:channelNum-1 {
					refname = mouseFile+"_$channel"+".ref"
					echo Loading reference file $refname
					ref = <;>
					iread ref refPath+mouseFile+"_$channel"+".ref"
					
					scalA = <;>
					echo
					echo Reading amplitude scalogram for $mouseFile channel $[channel+1]
					iread scalA scalAPath+mouseFile+'_$WTparamStr'+'$channel'+'.scalA'
					
					echo Analyzing harmonic amplitudes
					harmonicVals2 scalA ref hamp startXcorr endXcorr A
					delete scalA
					iwrite hamp hampPath+mouseFile+'_$RoSparamStr'+'$channel'+'.hamp'
					
					scalP = <;>
					echo
					echo Reading phase scalogram for $mouseFile channel $[channel+1]
					iread scalP scalPPath+mouseFile+'_$WTparamStr'+'$channel'+'.scalP'

					echo Analyzing phases of putative harmonics
					# Populates hphz with phase values at frequency harmonics
					harmonicVals2 scalP ref hphz startXcorr endXcorr P
					delete scalP
					iwrite hphz hphzPath+mouseFile+'_$RoSparamStr'+'$channel'+'.hphz'

					# CALCULATES N:M PHASE DIFFERENCES
					echo
					echo Calculating n:m phase differences				
					phaseDiff = Zero(hphz.ncol)
					foreach H 2:4 {
						phaseDiff = [pd_s hphz[0;:] hphz[H-1;:] H]
						nm = '1$H'
						echo Writing 1:$H phase difference signal
						write phaseDiff psyncPath+mouseFile+'_$RoSparamStr'+'$channel'+'.pd$nm'
					}
				echo DONE $mouseFile channel $[channel+1]
				echo ~~~~~~~~~~~~
				echo
			}

			echo DONE $mouseFile
			echo ***
			echo
			delete mouseFile
			delete hamp
			delete hphz

		}
	}

	if (PhaseDistHist==1) {
		echo 
		echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		echo Starting phase difference distribution generation
		source numdur2
		source pd3d_v3
		source addResult
		
		phase3Ddist = {}
		foreach i 0:359 {phase3Ddist+= Zero(360,360)}
# p0={} ;; sumpd=Zero(360) ;; mxpd=Zero(360) ;; foreach i 0:359 {p0+=<;> ;; iread p0[i] pd3dPath+"raw/"+"pd_matrix."+"$i" -a 360 360 ;; sumpd[i]=sum(p0[i]) ;; mxpd[i]=max(p0[i])} ;; foreach i 0:359 {p0[i]/=max(mxpd)}

# p1={} ;; sumpd=Zero(360) ;; mxpd=Zero(360) ;; foreach i 0:359 {p1+=<;> ;; iread p1[i] pd3dPath+pdDir+"/raw/"+"pd_matrix."+"$i" -a 360 360 ;; sumpd[i]=sum(p1[i]) ;; mxpd[i]=max(p1[i])} ;; foreach i 0:359 {p1[i]/=max(mxpd)}
		
		fNames = [file list lwrPath+"*.lwr"]
		fNames = [listv sort fNames]
		
		if (partialSWD==1) {
			pdDir = "partialSWDs_"+"$partialSWDrng[@<]"+"_"+"$partialSWDrng[@>]"
			if ([file list pd3dPath+pdDir].size==0) {
				file createdir pd3dPath "$pdDir"
			}
			pd3dPath = pd3dPath+pdDir+"/"
		}
		
		foreach session 0:fNames.size-1 {
			mouseFile = fNames[session]
			mouseFile[[str substr mouseFile '.lwr']] := ''
			
			annotationSig = <>
			read annotationSig annotPath+mouseFile+'.tannot'
			
			if (partialSWD==1) {
				annotRngs = [numdur2 [copy annotationSig] 1 swd]
				modAnnot = Zero(annotationSig.size)
				foreach i 0:annotRngs.size-1 {
					modStart = partialSWDrng[@<]+annotRngs[i][ranges][@<]
					modEnd = partialSWDrng[@>]+annotRngs[i][ranges][@<]
					if (modStart<annotRngs[i][ranges][@>] && modEnd<annotRngs[i][ranges][@>]) {
						modAnnot[modStart:modEnd]:=2
					}
				}
				annotationSig = modAnnot
				modRngs = [numdur2 [copy annotationSig] 1 swd]
				partStr = "$thisFolder -- $modRngs.size out of $annotRngs.size total SWDs used in $mouseFile"
				addResult drivePath partStr "$pdDir"+"_strain_info.txt"
			}
			
			sessionTotal = fNames.size
			thisSession = session+1
			
			echo
			echo Starting $mouseFile

			foreach channel channelRanges { 
				echo channel $channel in $mouseFile ($thisSession of $sessionTotal)

				pd12 = <>
				pd13 = <>
				pd14 = <>
				read pd12 psyncPath+mouseFile+'_'+'$RoSparamStr'+'$channel'+'.pd12'
				if (max(pd12)<=180 && min(pd12)<0) {pd12+=180}
				read pd13 psyncPath+mouseFile+'_'+'$RoSparamStr'+'$channel'+'.pd13'
				if (max(pd13)<=180 && min(pd13)<0) {pd13+=180}
				read pd14 psyncPath+mouseFile+'_'+'$RoSparamStr'+'$channel'+'.pd14'	
				if (max(pd14)<=180 && min(pd14)<0) {pd14+=180}

				phaseDiffs_x_time = Zero(4,pd12.size)
				
				phaseDiffs_x_time[0;:] = [copy pd12]
				delete pd12
				phaseDiffs_x_time[1;:] = [copy pd13]
				delete pd13
				phaseDiffs_x_time[2;:] = [copy pd14]
				delete pd14
				
				if (pd3dSWCs==1) {
					hphz = <;>
					iread hphz hphzPath+mouseFile+'_'+RoSparamStr+'$channel'+'.hphz'
					phaseDiffs_x_time[3;:] = der(hphz[0;:])*-1
					pd3dPath = pd3dPath+"pd3dSWC/"
					pd3d_v3 phase3Ddist phaseDiffs_x_time annotationSig -add_by_swc pd3dPath
				} else {
					pd3d_v3 phase3Ddist phaseDiffs_x_time annotationSig -add pd3dPath
				}
			}
			echo Phase difference distribution for $mouseFile DONE
			echo ~~~~~~~~~~~~~~~~~~
		}

		# To write phase degree matrices after all sessions complete
		# The -save_fiji option does not use phase difference x time &image nor any filter (only -add OPTION requires these arguments)
		# depending on whether LW running on windows or linux, pd3d action will have difference path formats, so to keep control at level of preproc append os onto pd3d action
		# (assuming it works) - alternately, just use if statements....
		if (os=='win') {
			pd3d_v3 phase3Ddist ~<0,0> <0> -save_fiji_win pd3dPath
		} elseif (os=='linux') {
			pd3d_v3 phase3Ddist ~<0,0> <0> -save_fiji_linux pd3dPath
		}
# p0 = {} ;; foreach i 0:359 {p0+=[copy phase3Ddist[i]]}
		echo DONE phase difference cube written to disk

		# to get ImageJ/Fiji output back into LW
	#	pd3d_v3 {} ~<0,0> <0> -load_fiji pd3dPath

	}

	if (PhaseDistSeq==1) {
		echo 
		echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		echo Starting write of phase difference trajectory sequences
		source numdur2
		source pd3d_v4
		source addResult
		source trajVel
		source freqInd
		
		# index to frequency
		ifc = [listv cread freqIndexFile F][0]
			
		fNames = [file list lwrPath+"*.lwr"]
		fNames = [listv sort fNames]

		foreach session 0:fNames.size-1 {
			mouseFile = fNames[session]
			mouseFile[[str substr mouseFile '.lwr']] := ''
			
			annotationSig = <>
			read annotationSig annotPath+mouseFile+'.tannot'
			annotRngs = [numdur2 [copy annotationSig] 1 swd]
			
			identifiers = {}
			identifiers.size = 3
			identifiers[0]=thisFolder
			identifiers[1]=mouseFile
			identifiers+=annotRngs
			
			echo 
			echo Starting $mouseFile - $annotRngs.size SWDs
	 
			foreach channel channelRanges {
				echo channel $channel in $mouseFile
				
				identifiers[2]=channel
				
				hphz = <;>
				phasefname = mouseFile+'_'+RoSparamStr+'$channel'+'.hphz'
				iread hphz hphzPath+phasefname
				fundFreq = der([copy hphz[0;:]])*-1
				delete hphz
				
				pd12 = <>
				pd13 = <>
				pd14 = <>
				read pd12 psyncPath+mouseFile+'_'+'$RoSparamStr'+'$channel'+'.pd12'
				read pd13 psyncPath+mouseFile+'_'+'$RoSparamStr'+'$channel'+'.pd13'
				read pd14 psyncPath+mouseFile+'_'+'$RoSparamStr'+'$channel'+'.pd14'	

				phaseDiffs_x_time = Zero(7,pd12.size)
				
				phaseDiffs_x_time[0;:] = [copy pd12]
				delete pd12
				phaseDiffs_x_time[1;:] = [copy pd13]
				delete pd13
				phaseDiffs_x_time[2;:] = [copy pd14]
				delete pd14

				# putting EEG into phaseDiffs_x_time for expediency
				rawEEG = <;>
				iread rawEEG lwrPath+mouseFile+'.lwr'
				phaseDiffs_x_time[3;:] = [copy rawEEG[channel;:]]
				delete rawEEG

				# putting indices for FF from ref file and actual FF from indices into phaseDiffs_x_time for expediency
				ref = <;>
				iread ref refPath+mouseFile+'_'+'$channel'+'.ref'
				phaseDiffs_x_time[4;:] = [copy ref[xcov_FF_indices;:]]
				phaseDiffs_x_time[5;:] = [freqInd [copy ref[xcov_FF_indices;:]] ifc]
				delete ref
				
				# putting power across fundamental frequency into phaseDiffs_x_time for expediency
				hamp = <;>
				iread hamp hampPath+mouseFile+'_'+RoSparamStr+'$channel'+'.hamp'
				phaseDiffs_x_time[6;:] = [copy hamp[0;:]]
				delete hamp
				
				pd3d_v4 identifiers phaseDiffs_x_time fundFreq -sequence pd3dPath
			}
			echo SWD trajectories written for $mouseFile DONE
			echo ~~~~~~~~~~~~~~~~~~
		}
	}

}

echo
echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
echo PROGRAM STARTED at $startt[0]:$startt[1]:$startt[2]
endt = [time current -l]
if (endt[1]<10) {endt[1] = '0$endt[1]'} else {endt[1] = '$endt[1]'}
if (endt[2]<10) {endt[2] = '0$endt[2]'} else {endt[2] = '$endt[2]'}
echo PROGRAM ENDED at $endt[0]:$endt[1]:$endt[2]
echo ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~