- MNE version: 1.0.2
- operating system: macOS Big Sur version 11.6.2
First of all I want to reproduce a result (a code) produce by normal python package, using men package.
So I would like to do a wavelet_morlet analysis to get the phase and then with the phase I can transform it to a complex version (on the unit cercle) by doing np.exp(1j*phase).
> epochsTFR_phase = tfr_morlet(epochs, freqs = freqs, n_cycles = n_cycles, return_itc=False, average=False,output='phase',use_fft = False,zero_mean=True)
> DA_epochsTFR_phase = np.transpose(epochsTFR_phase.data,(2,0,3,1)) ##(nTrials,nSensors,nfreq,nTimes) to (nfeq,nTrials,nTimes,nSensors)
> DA_epochsTFR_PHcomplex = np.exp(1j*DA_epochsTFR_phase)
Then I plot the cosine similarity (so the ratio between DA_epochsTFR_PHcomplex and my model that is also in complex)
I should have something like that
but I get something like that
The code that I want to reproduce, calculate by “hand” the wavelet like that
> def make_wavelet_CTS(self,f): #times,sensors for resting state data
> #one cycle at this frequency times by number of cycles in wavelet
> fi_range = int(self.p.N_cycles*1000.0/(self.p.freqs[f]*self.p.tdel))
> #half the wavlet cycles at lowest frequency minus half the wavelet cycles at this frequency
> PH_pad_at_f = self.p.PH_unpad_at_minf - int(self.p.N_cycles*500.0/(self.p.freqs[f]*self.p.tdel)) #zero padding at min frequency
> W = self.hamming(fi_range)
> f_indexes = []
> #how many sets of indexes to make? one for each sample in phase
> for s in range(self.p.PHSamples):
> #make a list of all data samples that are used in calculating each phase estimate
> #f_indexes belong to LTSSamples, so we move to the correct sample within TSfile, then pad
> padded_range = range(s+PH_pad_at_f,s+PH_pad_at_f+fi_range)
> f_indexes += padded_range
> streamCTS = self.source
> lpower = []
> phase = []
> for l in progressbar(range(self.p.nTrials), "Doing trial:", self.p.nTrials):
> lpower_sw = []
> phase_sw = []
> streamTS = streamCTS[l,:,:]
> for s in range(self.p.nSensors):
> streamTW = streamTS[f_indexes,s].reshape(self.p.PHSamples,fi_range) #tw
> base = streamTW.mean(1) #over w
> streamTW = streamTW - base[:,np.newaxis]
> t_f = 2.0*np.pi*self.p.N_cycles*np.arange(fi_range)/fi_range
> t_phi = np.exp(1j*(t_f)) * W
> f_com = streamTW * t_phi[np.newaxis,:]
> fourier = f_com.sum(1) #over w
> lpower_sw += [np.log10(np.absolute(fourier))[...,np.newaxis]]
> phase_sw += [np.angle(fourier)[...,np.newaxis]]
> lpower += [np.concatenate(lpower_sw,axis=-1)]
> phase += [np.concatenate(phase_sw,axis=-1)]
> return np.array(phase),np.array(lpower) #CTS
>
>
> def gaus(self,x,a,x0,sigma):
> return a*np.exp(-(x-x0)**2/(2*sigma**2))
> def hamming(self,x): #modified to multiple length window to calculate best fit
> from scipy.optimize import curve_fit
> if type(x) is int:
> n = x
> x = np.arange(n*5)
> else:
> n = len(x)*5
> print('hann window size',n)
> hann = 0.5*(1.0 - np.cos(2*np.pi*x/(n-1)))
> mean = sum(x*hann)/n
> sigma = np.sqrt(sum(hann*(x-mean)**2)/n)
> popt,pcov = curve_fit(self.gaus,x,hann,p0=[1,mean,sigma])
> G = self.gaus(x,*popt)
> return G[3::5] #reduce the resolution of the estimate
I want to have the same output results with men python, unfortunately when I plot the graphics it is far from being the same plots… (only the output of this functions make_wavelet_CTS and tfr_morlet differs)
So I don’t know how to get the same output with tfr_morlet arguments … Maybe the tfr_morlet is not constructed on the same calculation as above ?
Thank you very much for any type of informations !!