calcPhase.py 文件源码

def calcPhase():
    ipctx1 = dict()
    ipctx1["ip"] = np.zeros((1,2001))   
    dt = p.params["dt"] 

    randSamp = np.random.randint(0,len(AllAs),1000) 
    # Phases of three nuclei
    phasMus = dict()
    Flags = []
    Flags.append("SWA") 
    phasTAvsTA = []
    phasTIvsSTN = []
    phasTAvsSTN = []
    for i,samp in enumerate(randSamp):
        #temp = np.zeros((8,len(freqs),len(fftfreq)/10+1))
        A = AllAs[samp]
        B = AllBs[samp]
        Rates = psGA.calcRates(Flags,1.0,A,B,False,ipctx1,iptau=p.params["iptau"])



        # Filter the signals first at SWA frequency, or the phase difference between cortical input signal and TA,TI not calculated correctly due to multiple frequencies in the signal
        #B, A = sciSig.butter(2,np.array([0.0001,0.0005]),btype='low')
        B, A = sciSig.butter(2,np.array([0.00005]),btype='low')
        taFilt = sciSig.filtfilt(B, A, Rates['ta'])#, padlen=150)
        tiFilt = sciSig.filtfilt(B, A, Rates['ti'])#, padlen=150)
        stnFilt = sciSig.filtfilt(B, A, Rates['stn'])#, padlen=150)
        tG = np.arange(0,len(ipctx1["ip"][0])+dt,dt)

        fftfreq = np.fft.fftfreq(len(taFilt),d=dt)
        fftta = np.fft.rfft(taFilt-np.mean(taFilt))
        fftti = np.fft.rfft(tiFilt-np.mean(tiFilt))
        fftstn = np.fft.rfft(stnFilt-np.mean(stnFilt)) 

        fftip = np.fft.rfft(Rates['ipctx'] -np.mean(Rates['ipctx']))
        maxta = np.where(np.abs(fftta)==np.max(np.abs(fftta)))[0]
        maxti = np.where(np.abs(fftti)==np.max(np.abs(fftti)))[0]
        maxstn = np.where(np.abs(fftstn) == np.max(np.abs(fftstn)))[0]
        maxip = np.where(np.abs(fftip) == np.max(np.abs(fftip)))[0]

        print "maxta,maxti,maxstn,maxip",maxta,maxti,maxstn,maxip


        phasTAvsTA.append(np.angle(fftta[maxta]/fftta[maxta]))

        phasTIvsSTN.append(np.mean([np.angle(fftti[maxti]/fftstn[maxti]),np.angle(fftti[maxstn]/fftstn[maxstn])] ))
        phasTAvsSTN.append(np.mean([np.angle(fftta[maxta]/fftstn[maxta]),np.angle(fftta[maxstn]/fftstn[maxstn])] ))
    phasMus["ta_ta"] = phasTAvsTA
    phasMus["stn_ti"] = phasTIvsSTN
    phasMus["stn_ta"] = phasTAvsSTN

    pickle.dump(phasMus,open(path5+"Phases.pickle","w"))