Dear Julia peeps,
I am trying to relate brain activity to an amplitude modulated tone stimulus. As in, does the peak of brain activity occur in synchrony with the phase of the modulated tone? HypothesisTests.jl has this function for a Rayleigh Test: ( HypothesisTests.jl/circular.jl at master · JuliaStats/HypothesisTests.jl (github.com)) but I am not sure from the documentation how to use it. Below is the code I’m using to find instantaneous phase of a 2 Hz amplitude modulation wave and relate a series of mock peak latencies to it. I then take the phase as orientation for a series of unit vectors and find the sum of vectors for the length of the resultant vector (vstrength). I want to compare that resultant vector length between conditions and I believe I should be able to do that with the Rayleigh Test. Any help would be appreciated on how to use this function properly.
(I’m using this publication as reference: Auditory Cortex Phase Locking to Amplitude-Modulated Cochlear Implant Pulse Trains (physiology.org))
using DSP, FFTW, HypothesisTests # prep the data which will fluctuate peaklat1 = Int.(ceil.(rand(100)*100)) # list of fake latencies peaklat1[peaklat1 .> 50] = peaklat1[peaklat1 .> 50].+150 # create more variance peaklat2 = Int.(ceil.(rand(100)*100)) # list of fake latencies clickfreq = 2 # stimulus frequency sr = 1000 # sampling rate function getAMwave(clickfreq=2,sr=1000) t = 0:1/sr:1 # time course wave = sin.((2*π*clickfreq*t).+3π/2) # AM wave starting at the bottom return wave end # generate the AM wave used and get instantaneous phase wave = getAMwave(clickfreq,sr) # function above wavephase = angle.(hilbert(wave)) # Compute phase of wave # get the phase at each time point of signal peak amplitude orientation1 = wavephase[peaklat1] # break vector into components i and j to compute sums, averages, and then the resultant vector i1 = cos.(orientation1) j1 = sin.(orientation1) sumi1 = sum(i1) sumj1 = sum(j1) avgi1 = sumi1/length(i1) # average of component i avgj1 = sumj1/length(j1) # average of component j vstrength1 = sqrt((avgi1^2)+(avgj1^2)) # find the vector based on i and j components (Pythagorean) # now the same for the second latency series orientation2 = wavephase[peaklat2] i2 = cos.(orientation2) j2 = sin.(orientation2) sumi2 = sum(i2) sumj2 = sum(j2) avgi2 = sumi2/length(i2) avgj2 = sumj2/length(j2) vstrength2 = sqrt((avgi2^2)+(avgj2^2)) # what I feel I should be able to do: RayleighTest(vstrength1,vstrength2) # or RayleighTest([vstrength1,vstrength2],100) # 100 observations per condition