AkvoFluo/GUI/quad_corr.py

# import matplotlib as mpl
# mpl.use("pgf")
# pgf_with_pdflatex = {
#     "pgf.texsystem": "pdflatex",
#     "pgf.preamble": [
#          r"\usepackage{amsmath}",
#          r"\usepackage{amssymb}",
#          r"\usepackage{timet}",
#          ]
# }
#          #r"\usepackage[utf8x]{inputenc}",
#          #r"\usepackage[T1]{fontenc}",
#          #r"\usepackage{cmbright}",
# mpl.rcParams.update(pgf_with_pdflatex)
# 
# from matplotlib import rc
# rc('font',**{'family':'serif','serif':['times']})
# ## for Palatino and other serif fonts use:
# #rc('font',**{'family':'serif','serif':['Palatino']})
# #rc('text', usetex=True)

#import matplotlib
#from matplotlib import rc

#matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{timet,amsmath}"]
#rc('font',**{'family':'serif','serif':['timet']})
#rc('text', usetex=True)
from GJIPlot import *

import numpy as np
import matplotlib.pyplot as plt
from smooth import smooth
#from invertColours import *
from decay import *
from scipy import signal
from notch import * 
import rotate # why is it called loss?

BLACK = False
labs = 16

if BLACK:
    labcol = 'white'
else:
    labcol = 'black'

def Kw(a, phi2, s, rT2, wL): 
    return  a*((np.sin(phi2)*s + ((1/rT2)*np.sin(phi2)) + wL*np.cos(phi2) ) / (wL**2+(s+1/rT2)**2 ))


def quadDetect(T, vL, wL, dt, xn, DT):
    # decimate
    # blind decimation
    # 1 instead of T 
    irsamp = (T) * int(  (1./vL) / dt) # real 
    iisamp =       int(  ((1./vL)/ dt) * ( .5*np.pi / (2.*np.pi) ) ) # imaginary
   

    dsamp =  int( DT / dt) # real 

    iisamp += dsamp

    ############################################################
    # simple quadrature-detection via sampling
    xr = xn[dsamp::irsamp]
    xi = xn[iisamp::irsamp]
    phase = np.angle( xr + 1j*xi )
    abse = np.abs( xr + 1j*xi )

    # times
    ta = np.arange(0, TT, dt)
    te = np.arange(DT, TT, TT/len(abse))

    #############################################################
    # hilbert transform
    ht = signal.hilbert(xn) #, 100))
    he = np.abs(ht)         #, 100))
    he = ht*np.exp( -1j*wL*t )  # ht*e^{-jwLt}  # phase envelope

    #hp = ((np.angle(ht[dsamp::irsamp]))) 
    #[Q, I, tt] = rotate.quadrature(T, vL, wL, dt, xn, DT, t)
    [E0,df,phi,T2] = quadratureDetect(np.imag(he), np.real(he), t)
    D = rotate.RotateAmplitude(np.real(he), np.imag(he), phi, df, t)
    he = D.imag 
    #print ("df", df, "phi", phi)
    """
    plt.plot(xn)
    plt.plot(np.real(he))
    plt.plot(np.imag(he))
    plt.plot(D.imag)
    plt.plot(D.real)
    print ("phi", phi)
    plt.plot( np.unwrap(np.arctan(np.unwrap(np.angle( he ) + np.pi/2 ))) )
#    plt.plot(np.imag(ht))
#    plt.plot(xn)

    # plot FFT

    XN  = np.fft.rfft(xn)
    plt.figure()
    plt.plot(XN.real, label='orig')
    #XNN = XN * np.exp(1j*phi-1.)
    """

    """
    phi = np.arange(-np.pi, np.pi, .001)
    best = 1e12
    bestPhi = 0
    for p in phi:
        XNN = XN * np.exp(1j*p)
        if ( (np.linalg.norm( XNN.real[XNN.real < 0] )) < best):
            best = np.linalg.norm( XNN.real[XNN.real < 0])
            bestPhi = p
        #print (len( XNN.real[XNN.real < 0] ))
        #best = min(best, -( np.sum(XNN.real[XNN.real < 0]) ))
        #print (XNN[XNN.real < 0]) 
        #plt.plot(XNN.real[XNN.real<0])
    print ("best", best, bestPhi)
    XNN = XN * np.exp(1j*bestPhi)
    """
   
    """ 
    plt.plot(XNN.real, label='phased')
    plt.plot(XNN.imag, label='imag')
    plt.legend()
    plt.figure()
    plt.plot(np.angle(XN))
    plt.plot(np.angle(XNN))
    plt.show()
    exit()
    """

    #############################################################
    # Resample ht, no gate integrate
    """
    htd = signal.decimate(he, 100, ftype='fir') 
    td = signal.decimate(t, 100, ftype='fir')
    #[htd, td] = signal.resample(he, 21, t) 
    #toss first, and use every third 
    htd = htd[1::5]
    td = td[1::5]
    """

    #############################################################
    # Resample ht and gate integrate to 20 gates per decate
    GPD = 20 # Gates per decade
    GI = 15  # initial 15 times are just used

    # first just resample to 20 Gates per decade
    #nrs = (np.log(t[-1]) - np.log10(t[0]))
    tdd = np.logspace( np.log10(16*dt), np.log10(t[-1]), 80-GI, base=10, endpoint=True) 
    tdl = tdd[0:-1]    # these are left edges
    tdr = tdd[1::]    # these are left edges
    td = (tdl+tdr) / 2.
    
    # Just use the first 20 datapoints
    td = np.concatenate( ([t[0:GI] ,  td]) )

    htd = np.zeros( len(td) )
    #htd[0:15] = he[0:15]
    htd[0:GI] = he[0:GI]

    ii = 15
    isum = np.zeros( len(td) )
    Vars = np.ones( len(td) ) * .15**2
    isum [0:15] = 1
    #Vars [0:15] = .15**2  # variance 
    #var = .15**2
    for itd in range(GI, len(he)):
        #print (itd)   
        if ( t[itd] > tdr[ii-GI] ):
            ii += 1
        isum[ii] += 1
        htd[ii] += he[ itd ]
        Vars[ii] += .15**2
    htd /= isum
    #Vars = ( Vars) / isum
    Stds = np.sqrt( Vars) / isum
    Vars = Stds**2 
    #plt.figure()
    #plt.gca().errorbar(  td, htd, fmt='o', yerr=Vars )
    #plt.gca().set_xscale('log', basex=10)
    
    #    htd = he[itd] # (int)(dt*itd) ]
    #td = []
    # OK, so tdd are the desired time gates, but we might not be able to satisfy all of them

        

    #htdt = signal.decimate(he, 10, ftype='fir') 
    ##tdt = signal.decimate(t, 10, ftype='fir')
    #tdt = t[::10]
    #[htd, td] = signal.resample(he, 21, t) 
    
    ##toss first, and use every third 
    #htd = [] # htdt[1::5]
    #td = tdt[1::15]
    #for ix in np.arange(1,len(htdt),15):
    #    if ix>15 and ix<len(htdt)-15:
    #        htd.append( np.mean(htdt[ix-7:ix+8]) ) # stupid could do recursive but who gives a shit
    #    else:
    #        htd.append( htdt[ix] ) # stupid could do recursive but who gives a shit
    
    ###############################################
    # Plots
    #plt.plot(ta, xn)
    #plt.plot(te, abse, '-.', linewidth=2, markersize=10)
    #plt.plot(ta, he, '.', markersize=10 )
    #plt.plot(td, htd, color='green', linewidth=2)
    # Phase Plots
    #ax2 = plt.twinx()
    #ax2.plot(te, hp, '.', markersize=10, color='green' )
    #ax2.plot(te, phase, '-.', linewidth=2, markersize=10, color='green')
    #plt.show()
    #exit()

    #####################################################################
    # regress raw signal
    
    #[peaks, times, ind] = peakPicker(xn, wL, dt)
    #[a0,b0,rt20] =  regressCurve(peaks,  times) #,sigma2=1,intercept=True):
    #plt.figure()
    #plt.plot(Vars)
    #plt.show()
    #exit()

    
    dsamp = 1 # int( DT / dt) # real 
    # regress analytic signal
    [b0,rt20] =  regressCurve(he,  t, intercept=False) #, sigma2=(1./Stds)/np.max(1/Stds)) #,sigma2=1,intercept=True):
    #[a0,b0,rt20] =  regressCurve(he,  t, intercept=True, sigma2=(1./Vars)/np.max(1/Vars)) #,sigma2=1,intercept=True):
    #[b0,rt20] =  regressCurve(he[dsamp::],  t[dsamp::], intercept=False) #,sigma2=1,intercept=True):
    #[a0,b0,rt20] =  regressCurve(he,  t) #,sigma2=1,intercept=True):
   
    # regress downsampled 
    #[a,b,rt2] =  regressCurve(abse,  t[dsamp::irsamp], intercept=True) #,sigma2=1,intercept=True):
    [b,rt2] =  regressCurve(htd,  td, intercept=False,  sigma2=(1./Stds)/np.max(1/Stds) ) #,sigma2=1,intercept=True):
    
    return irsamp, iisamp, htd, b0, rt20, ta, b, rt2, phi, td, he, Stds
    #return irsamp, iisamp, abse, a0, b0, rt20, times, a, b, rt2, phase

if __name__ == "__main__":
    
    dt = 1e-4
    TT = 1.5
    t = np.arange(0, TT, dt)
    vL = 2000.
    wL = 2.*np.pi*vL
    wL2 = 2.*np.pi*(vL+10 ) # 3 Hz off 
    zeta = -np.pi/2 # .0 #np.pi/2.21
    t2 = .050

    # TODO, show boxplots of error of regressions parameters, and randomly sample T2 and Vn0 for realistic parameters. 
    # do a lot of random realizations.   

    N = 5000
    T2 =  [.03, .075, .125, .2, .250] # <=== Final 
    #T2 =  [.3]
    NT2 = len(T2)
    RALL = np.zeros((N, NT2))
    RDEC = np.zeros((N, NT2))
    BALL = np.zeros((N, NT2))
    BDEC = np.zeros((N, NT2))
    RF = np.zeros((N, NT2))
    AF = np.zeros((N, NT2))
    RM = np.zeros((N, NT2))
    AM = np.zeros((N, NT2))
    RFC = np.zeros((N, NT2))
    AFC = np.zeros((N, NT2))
    
    sT2 = [str(int(1e3*i)) for i in T2]
    it2 = 0 
    for t2 in T2:
        print(("simulating", t2))


        xe = np.exp(-t/t2)    
 
        for it in range(N):
        
            zeta = (np.random.rand() - .5) #np.random.rand()
            #print ("zeta", zeta)
            # TODO, right here do distribution of phases  
            xs = np.exp(-t/t2) * np.sin(wL*t + numpy.pi*zeta)  #\
              # + np.exp(-t/t2) * np.sin(wL*t  ) 
            
            # pretty low, highly-correlated noise
            #xn =  ( (smooth(np.random.normal(0,.4,len(xs)) + \
            #    np.random.random_integers(-1,1,len(xs))*0.5*np.random.lognormal(0, .25, len(xs)) + \
            #    np.random.random_integers(-1,1,len(xs))*.0002*np.random.weibull(.25, len(xs)), 60)))
            #xn += np.random.normal( 0, .02,len(xs) ) 
            #print(np.std(xn))

            xc = np.random.normal(0,.01,len(xs))  + (np.sign(xs) *  (smooth(np.random.normal(0,.25,len(xs)+160) + \
                np.random.random_integers(-1,1,len(xs)+160)*0.15*np.random.lognormal(0, .25, len(xs)+160) + \
                np.random.random_integers(-1,1,len(xs)+160)*.0002*np.random.weibull(.25, len(xs)+160), 80))[80:-80] )           
            
            #print(np.std(xc))
            
            xn = xs + xc 


            #xn += xs            

            #nn = np.random.normal( 0, .15,len(xs) ) 
            #xn = xs + nn # np.random.normal(0,.15,len(xs)) 
            
            #xn = xs + (np.sign(xs) *  np.absolute(smooth(np.random.normal(0,.5,len(xs))) )) #+ \
                #np.random.random_integers(-1,1,len(xs))*0.6*np.random.lognormal(0, .35, len(xs)) + \
                #np.random.random_integers(-1,1,len(xs))*.004*np.random.weibull(.25, len(xs)), 60)))


            ###################################################################
            # Bandpass filter, concatenate the two and run both directions
            #xraw = xn
            #plt.plot(xn)
            
            #xf = notch( xn[::-1], .925, 2000., dt )            
            #xn = (xn[::-1]-xf)[::-1]
            
            #xf = notch( xn[::-1], .9, 2000., dt )            
            #xn = (xn[::-1]-xf)[::-1]

            """ 
            print("sigma xf", np.std(xf))
            print("sigma nn", np.std(nn))
            print("sigma xn", np.std(xn[-100::]))
            #plt.plot(xraw)
            plt.plot(xn)
            plt.plot(xs)
            #plt.plot(xf)
            #plt.plot(xn[::-1]-xf)
            plt.show()
            exit()
            """

            # quadrature detection downsampling
            T = 50    # sampling period, grab every T'th oscilation
            DT = 0. #.002 #85  # dead time ms
            [irsamp, iisamp, abse, b0, rt20, times, b, rt2, phase, tdec, he, Stds] = quadDetect(T, vL, wL, dt, xn, DT)
            RALL[it, it2] = 1e2 * ((rt20-t2)/t2)
            BALL[it, it2] = 1e2*(b0 -1.)
            BDEC[it, it2] = 1e2*(b - 1.)
            RDEC[it, it2] = 1e2 * ((rt2-t2)/t2)

            # Fourier domain fitting
            dsamp = int( DT / dt) # real 
            #X = np.fft.rfft(xn[dsamp::]) #* dt
            #V = np.fft.fftfreq(len(xn[dsamp::]), dt)
            X = np.fft.rfft(xn) #* dt
            V = np.fft.fftfreq(len(xn), dt)

            #V[len(X)-1] *= -1.
            w = V * 2.*np.pi
 
            # Just do a window
            nf = 50 # 150 =  200 Hz window
            ws = np.abs(w[0:len(X)])
            vs = ws/(2.*np.pi)
            iL =  vL / vs[1]

            # Do complex fit 
            [af, bf] = regressModulus(ws[iL-nf:iL+nf] , wL, dt*X[iL-nf:iL+nf] )
            AM[it, it2] =  1e2*(af - 1.) # backwards extrapolate
            RM[it, it2] = 1e2 * ((bf-t2)/t2)

            try:
                [af, bf, nb] = regressSpecComplex(ws[iL-nf:iL+nf], wL, dt*X[iL-nf:iL+nf], known=False)
            except:
                print("munged FD regression")
                af = 0
                bf=10.
                nb = 0
    
            cr = Kw(af, nb, 1j*ws[iL-nf:iL+nf], bf, wL) 
            
            AFC[it, it2] =  1e2*(af - 1.) # backwards extrapolate
            afe = af #  (af * np.exp(DT/bf)) #, af # backwards extrapolate
            RFC[it, it2] = 1e2 * ((bf-t2)/t2)
            """
            # rotate phase
            phi = np.arange(0, 2.*np.pi, .001)
            best = 1e12
            bestPhi = 0
            for p in phi:
                XNN = X * np.exp(1j*p)
                if ( (np.linalg.norm( XNN.real[XNN.real < 0] )) < best):
                    best = np.linalg.norm( XNN.real[XNN.real < 0])
                    bestPhi = p
            nb2 = bestPhi
            nb3 = nb
            """
            #print ("nb", nb, "zeta", zeta*pi, "nb2", nb2, np.pi/2-zeta*pi)
            Xc = X[iL-nf:iL+nf]
            
            #X = (X * np.exp(1j*(np.pi/2-zeta*pi) )).real # Correct
            X = (X * np.exp(1j*(np.pi/2-nb) )).real
            X = X[iL-nf:iL+nf]
            #print (phase, bestPhi, phase-bestPhi) 
            #plt.plot(ws, np.abs(X))
            #plt.show()
            #exit()           

            #[a, b] = regressSpec(w[0:len(X)], wL, dt*X)
            #[af, bf, nb] = regressSpecComplex(ws[iL-nf:iL+nf], wL, dt*X[iL-nf:iL+nf])
            #[af, bf, nb] = regressSpec(ws[iL-nf:iL+nf], wL, dt*X[iL-nf:iL+nf])
            [af, bf, nb] = regressSpec(ws[iL-nf:iL+nf], wL, dt*X)
            #AF[it, it2] =  af  - 1.
            #AF[it, it2] =  1e2*(af * np.exp(DT/bf)  - 1.) # backwards extrapolate
            AF[it, it2] =  1e2*(af - 1.) # backwards extrapolate
            afe = af #  (af * np.exp(DT/bf)) #, af # backwards extrapolate
            RF[it, it2] = 1e2 * ((bf-t2)/t2)
            
#             print af, bf, nb
# 
#             plt.figure()
#             plt.plot(xn)           
# 
#             plt.figure() 
#             plt.plot(V[0:len(X)], dt * np.abs(X), label='mod')
#             xp = af* np.abs( wL / (wL**2 + (-1j*ws[iL-nf:iL+nf]  + 1./bf)**2 ) )
#             xpc =  np.abs( wL / (wL**2 + (-1j*ws[iL-nf:iL+nf]  + 1./t2)**2 ) )
#             #plt.plot(vs[iL-50:iL+50], dt*np.real(X[iL-50:iL+50]))
#             #plt.plot(vs[iL-50:iL+50], dt*np.imag(X[iL-50:iL+50]))
#             plt.plot(vs[iL-nf:iL+nf], xp, label='predicted')
#             plt.plot(vs[iL-nf:iL+nf], xpc, label='actual')
#             plt.legend()

        it2 += 1

    # a little plot
    if BLACK:
        quadfig = plt.figure(0, facecolor='black', figsize=[pc2in(20), pc2in(32)])
    else: 
        quadfig = plt.figure(0, figsize=[ pc2in(20), pc2in(32) ])
    qa = quadfig.add_axes([.1, .6, .8,  .375])

    if BLACK:
        qa.plot(1e3*t, xn, label = 'noisy signal', color='yellow')
        invertColours(qa)
    else:
        qa.plot(1e3*t, xn, color='blue', alpha=.25)
        #qa.plot(t, xn, label = 'noisy signal', color='blue')

    #fa.plot(X, V)


    #plt.plot(t[::irsamp], xn[::irsamp], '.-', marker = 'b', markersize=18, label='real envelope')
    #plt.plot(t[iisamp::irsamp], xn[iisamp::irsamp], '.-', marker = 'imaginary envelope', markersize=18, label='imaginary envelope')
    deSpine(plt.gca())

    if BLACK:
        l2 = plt.plot(1e3*times, he, color='red', label='CA')
        l1 = plt.plot(1e3*tdec, abse, '.', color='green', markersize=12, label='time gates')
        leg = plt.legend( frameon=False, scatterpoints=1, labelspacing=0.2, numpoints=1 )
        invertLegend(leg)
    else:
        l2 = qa.plot(1e3*times, he, color='red') #, label='CA')
        #l1 = qa.plot(1e3*tdec, abse, 'o', color='green', markersize=4, label='time gates')

        #leg = plt.legend( frameon=True, scatterpoints=1, numpoints=1, labelspacing=0.25 )
        #rect = leg.get_frame()
        #rect.set_facecolor(light_grey)
        #rect.set_linewidth(0.0)
        #rect.set_alpha(0.5)

    #qa.set_xlabel("time [s]")
    #qa.set_title("time-domain")
    
    # regress curve using all times
    re0 = b0 * np.exp(-times/rt20)

    # regress curve based on decimated data
    irsamp = 1
    re = b * np.exp(-t[::irsamp]/rt2)   
        
    qa.plot(1e3*t, xe, label = 'true envelope', linewidth=2, color='blue')
    #qa.plot(1e3*t[::irsamp],re, linewidth=2, label='gated regression', color='green')
    qa.plot(1e3*t[::irsamp],re, linewidth=1, color='green')
    qa.plot(1e3*times,re0, '--',linewidth=1, label='CA regression', color='black')
    l1 = qa.errorbar( 1e3*tdec, abse, yerr=Stds, fmt='o', color='green', markersize=3, label='time gates (1 $\sigma$ EB)')
    qa.set_xscale('log', basex=10)
        
 
    leg = qa.legend(frameon=True, scatterpoints=1, numpoints=1, labelspacing=0.2)
    fixLeg(leg)


#       Frequency-domain axis, did not bring much information. So removed. 
    #ff = plt.figure()
    #fa = quadfig.add_axes([.1, .05,  .8, .25])
    fa = quadfig.add_axes([.1, .1, .8, .375])
 
    #ax2 = plt.gca().twinx()
    #l2 = ax2.plot(t[iisamp::irsamp], phase, '.-', marker = 'imaginary envelope', markersize=18, label='phase')
    #plt.legend([l1,l2])
    #invertColours(qa)
    
    #fa.plot(V[0:len(X)], dt * np.abs(X), label='mod')

    #fmla = robjects.Formula('XX ~ a*(.5/rT2)  / ((1./rT2)^2 + (w-wL)^2 )')
    #xp = af * np.abs( wL / (wL**2 + (-1j*ws[iL-nf:iL+nf]  + 1./bf)**2 ) )
    #xpc =     np.abs( wL / (wL**2 + (-1j*ws[iL-nf:iL+nf]  + 1./t2)**2 ) )



    xp = af * ( (.5/bf) / (  ((1./bf)**2 + (ws[iL-nf:iL+nf]-wL)**2 ) ) )
    xpc =     ( (.5/t2) / (  ((1./t2)**2 + (ws[iL-nf:iL+nf]-wL)**2 ) ) )
    #print ("af", af) 
    #plt.plot(vs[iL-50:iL+50], dt*np.real(X[iL-50:iL+50]))
    #plt.plot(vs[iL-50:iL+50], dt*np.imag(X[iL-50:iL+50]))
    #fa.plot(np.abs(V[0:len(X)]), dt * np.abs(X), color='red', linewidth=2, label='FT modulus')
    fa.plot(vs[iL-nf:iL+nf], dt*X, color='red', linewidth=2, label='phased')
    fa.plot(vs[iL-nf:iL+nf], xp, '--',linewidth=1, label='regressions', color='black')
    #fa.plot(vs[iL-nf:iL+nf], dt*Xc.real, color='red', '--', linewidth=2, label='real')
    fa.plot(vs[iL-nf:iL+nf], dt*Xc.real, '-', color='blue', linewidth=2, label=r'$\mathrm{Re}$ FFT ')
    fa.plot(vs[iL-nf:iL+nf], dt*Xc.imag, '--', color='blue', dashes=(1.5,.75), linewidth=2, label=r'$\mathrm{Im}$ FFT')
    
    fa.plot(vs[iL-nf:iL+nf], cr.real, '--', color='black', linewidth=1) #, label=r'reg')
    fa.plot(vs[iL-nf:iL+nf], cr.imag, '--', color='black', dashes=(1.5,.75), linewidth=1) #, label=r'reg')
    #ax.plot(dead_freqs[0:ni], np.imag(WSINDP2[0:ni]), "--", dashes=(1.5,.75),  color='blue', linewidth=1, label=r"$\mathrm{Im} ( \mathrm{DP} )$")
    #fa.plot(vs[iL-nf:iL+nf], xpc, linewidth=2, label='actual', color='blue')

    leg = fa.legend(frameon=True, scatterpoints=1, numpoints=1, labelspacing=0.2)
    fixLeg(leg)

    fa.set_xlabel("frequency [Hz]")
    #fa.set_title("frequency-domain")
    fa.set_xlim([1950,2050])
    #leg = plt.legend()

    #rect = leg.get_frame()
    #rect.set_facecolor(light_grey)
    #rect.set_linewidth(0.0)
    #rect.set_alpha(0.5)
    
    #ra = quadfig.add_axes([.1, .1,  .8, .25])
    #vs = .0125
    #ra = quadfig.add_axes([.1, .1-vs, .8, .375])
    #ra.set_title("decay envelope regressions")

    #plt.savefig('quadri.png',dpi=200, facecolor='black', edgecolor='black')

    # real imaginary plot
    #rfig = plt.figure(facecolor='black')
    #ra = rfig.add_axes([.1,.1,.8,.8])
    #ra.plot(t, xn, label = 'noisy signal', color='yellow')
    #ra.plot(t[::irsamp], xn[::irsamp], '.-', marker = 's', markersize=10, label='real signal')
    #ra.plot(t[iisamp::irsamp], xn[iisamp::irsamp], '.-', marker = 'v', markersize=10, label='imaginary envelope')
    #ra.plot(tdec, abse, '.-', color='green', markersize=18, label='decimated modulus')
    #ra.plot(t[iisamp::irsamp], phase, '.-', color='green', markersize=18, label='decimated phase')
    #leg = plt.legend()
    #invertLegend(leg)
    #invertColours(ra)
    #plt.savefig("quad2.png", dpi=200, facecolor='black', edgecolor='black')



    plt.figure(0) 
    deSpine(plt.gca())
    if BLACK:
        ra.plot(t[iisamp::irsamp],re, linewidth=2, label='decimated regression', color='green')
        ra.plot(times,re0, linewidth=2, label='all data regression', color='cyan')
        ra.plot(t, xe, label = 'true', linewidth=4, color='magenta')
        ra.set_xlabel("time", fontsize=16, color='white')
        leg = plt.legend( labelspacing=0.2, scatterpoints=1,frameon=False )
        fixLeg(leg)
        invertLegend(leg)
        plt.savefig('noisytime.pdf',dpi=600, facecolor='black', edgecolor='black')
    else:

        #ra.plot(1e3*times, afe*np.exp(-times/bf), '-', color='purple', linewidth=2, label='FD modulus')
        qa.set_xlabel("time [ms]") #, fontsize=16)
        

        #rect = leg.get_frame()
        #rect.set_facecolor(light_grey)
        #rect.set_linewidth(0.0)
        #rect.set_alpha(0.5)
        qa.set_ylim([-1.25,3.25])
        plt.savefig('noisytime_corr.pdf',dpi=600)
        plt.savefig('noisytime_corr.eps',dpi=600)

    #plt.show()
    #exit()
    
    #plt.savefig('quadreg.png',dpi=200, facecolor='black', edgecolor='black')

    # Just show envelopes
    #efig = plt.figure(facecolor='black')
    #ea = efig.add_axes([.1,.1,.8,.8])
    #ea.plot(t[iisamp::irsamp],re, linewidth=2, label='decimated regression', color='yellow')
    #ea.plot(times,re0, linewidth=2, label='all data regression', color='cyan')
    #ea.plot(t, xe, label = 'clean envelope', linewidth=2, color='magenta')
    #leg = plt.legend()
    #invertLegend(leg)
    #invertColours(ea)
    #plt.savefig("envelopes.png", dpi=200, edgecolor='black', facecolor='black')    

    # boxplot
    if BLACK:
        boxfig = plt.figure( facecolor='black', edgecolor='black')
    else:
        boxfig = plt.figure( figsize=[pc2in(20), pc2in(32)] )
    ##############
    #  b2    r2  #
    #  b1    r1  #
    #  b0    r0  #
    ##############
    vs =  0 #-3*.0125
    vs2 = 0 #-2*.0125

    # Time domain    
    b1 = boxfig.add_axes([.15,.075, .35, .15]) # bottom we want labels 
    r1 = boxfig.add_axes([ .6,.075, .35, .15])
    
    b0 = boxfig.add_axes([.15,.25, .35, .15])
    r0 = boxfig.add_axes([ .6,.25, .35, .15])
    
    # Frequency domain
    b2 = boxfig.add_axes([.15,.425, .35, .15])
    r2 = boxfig.add_axes([ .6,.425, .35, .15])
    
    b4 = boxfig.add_axes([.15,.6, .35, .15])
    r4 = boxfig.add_axes([ .6,.6, .35, .15])
    
    b3 = boxfig.add_axes([.15,.775,  .35, .15])
    r3 = boxfig.add_axes([ .6,.775,  .35, .15])
    
    for ax in [b1,r1,b0,r0,b2,r2,b3,r3,b4,r4]:
        deSpine(ax)
    
    B0 = b0.boxplot( BALL , sym='', bootstrap=5000)# , notch=True )
    plt.setp(b0, 'xticklabels', [])
    
    R0 = r0.boxplot( RALL , sym='', bootstrap=5000)# , notch=True )
    plt.setp(r0, 'xticklabels', []) 
    
    B1 = b1.boxplot( BDEC , sym='', bootstrap=5000)# , notch=True )
    plt.setp(b1, 'xticklabels', sT2) 
    
    R1 = r1.boxplot( RDEC , sym='', bootstrap=5000)# , notch=True )
    plt.setp(r1, 'xticklabels', sT2) 
    
    B2 = b2.boxplot( AF , sym='', bootstrap=5000)# , notch=True )
    R2 = r2.boxplot( RF , sym='', bootstrap=5000)# , notch=True )
    plt.setp(r2, 'xticklabels', []) 
    plt.setp(b2, 'xticklabels', []) 
    
    B3 = b3.boxplot( AFC , sym='', bootstrap=5000)# , notch=True )
    R3 = r3.boxplot( RFC , sym='', bootstrap=5000)# , notch=True )
    plt.setp(r3, 'xticklabels', []) 
    plt.setp(b3, 'xticklabels', []) 
    
    B4 = b4.boxplot( AM , sym='', bootstrap=5000)# , notch=True )
    R4 = r4.boxplot( RM , sym='', bootstrap=5000)# , notch=True )
    plt.setp(r4, 'xticklabels', []) 
    plt.setp(b4, 'xticklabels', []) 
    
    b1.set_xlabel(r"true $T_2^*$ [ms]", color=labcol) #, fontsize=12) 
    r1.set_xlabel(r"true $T_2^*$ [ms]", color=labcol) #, fontsize=12) 
    
    b2.set_ylabel(r"phased FD", color=labcol) #, fontsize=12) 
    b3.set_ylabel(r"complex FD", color=labcol) #, fontsize=12) 
    b4.set_ylabel(r"modulus FD", color=labcol) #, fontsize=12) 
    b0.set_ylabel(r"CA", color=labcol) #, fontsize=12) 
    b1.set_ylabel(r"gated CA", color=labcol) #, fontsize=12) 
    
    #r0.set_ylabel(r"parameter estimate error") 
    #b1.set_xlabel(r"$T_2^*$ decay parameter [ms]") 

    if BLACK:
        for box in [B0, R0, B1, R1]: 
            box['boxes'][0].set_color('w')
            box['boxes'][0].set_linewidth(2)
            box['medians'][0].set_color('y')
            box['medians'][0].set_linewidth(2)
            box['fliers'][0].set_color('w')
            box['fliers'][1].set_color('w')
            box['whiskers'][0].set_color('w')
            box['whiskers'][1].set_color('w')
#     else:
    #it2 = 0
    #for ax in B0, r0, b1, r1: 
        #for line in ax.xaxis.get_ticklabels():
        #        print  line.get_text()
        #        line.set_text(str(T2[it2]))
        #        print  line.get_text()
        #setp(gca(), 'xticklabels', ['first', 'second', 'third']) 
                #line.set_fontsize(16)   
     #   it2 += 1 
#             box['boxes'][0].set_color('red')
#             box['boxes'][0].set_linewidth(3)
#             box['medians'][0].set_color('y')
#             box['medians'][0].set_linewidth(3)
#             box['fliers'][0].set_color('w')
#             box['fliers'][1].set_color('w')
#             box['whiskers'][0].set_color('w')
#             box['whiskers'][1].set_color('w')
        

    #r0.set_title(r"analytic signal $T_2^*$", color=labcol)
    #b0.set_title(r"analytic signal $V_N^{(0)}$", color=labcol)
    r3.set_title(r"$\tilde{T}_2^*$ relative error [\%]", color=labcol) #, fontsize=12)
    b3.set_title(r"$\tilde{V}_N^{(0)}$ absolute error [\%]", color=labcol) #, fontsize=12)

    if BLACK:
        invertColours(b0)
        invertColours(r0)
        invertColours(b1)
        invertColours(r1)
        #plt.savefig("boxplot.png", dpi=200, facecolor='black', edgecolor='black')
    #else:
    #plt.figure()
    #plt.boxplot( BDEC, bootstrap=1000 )
    #plt.title("dec")    
        #plt.savefig("boxplot.png", dpi=200)

    bylow = min(b0.get_ylim()[0], b1.get_ylim()[0])
    bylow = min(bylow, b2.get_ylim()[0])
    bylow = min(bylow, b3.get_ylim()[0])
    bylow = min(bylow, b4.get_ylim()[0])

    byhi = max(b0.get_ylim()[1], b1.get_ylim()[1])
    byhi = max(byhi, b2.get_ylim()[1])
    byhi = max(byhi, b3.get_ylim()[1])
    byhi = max(byhi, b4.get_ylim()[1])
    
    rylow = min(r0.get_ylim()[0], r1.get_ylim()[0])
    rylow = min(rylow, r2.get_ylim()[0])
    rylow = min(rylow, r3.get_ylim()[0])
    rylow = min(rylow, r4.get_ylim()[0])

    ryhi = max(r0.get_ylim()[1], r1.get_ylim()[1])
    ryhi = max(ryhi, r2.get_ylim()[1])
    ryhi = max(ryhi, r3.get_ylim()[1])
    ryhi = max(ryhi, r4.get_ylim()[1])
 
    b0.set_ylim((bylow, byhi))
    b1.set_ylim((bylow, byhi))
    b2.set_ylim((bylow, byhi))
    b3.set_ylim((bylow, byhi))
    b4.set_ylim((bylow, byhi))
    
    r0.set_ylim((rylow, ryhi))
    r1.set_ylim((rylow, ryhi))
    r2.set_ylim((rylow, ryhi))
    r3.set_ylim((rylow, ryhi))
    r4.set_ylim((rylow, ryhi))

    
    if BLACK:
        #plt.savefig("boxplotlim2.pdf", dpi=600, facecolor='black', edgecolor='black')
        plt.savefig("boxplotlim_2.pdf", dpi=600, facecolor='black', edgecolor='black')
        plt.savefig("boxplotlim_2.eps", dpi=600, facecolor='black', edgecolor='black')
    else:
        #plt.savefig("boxplotlim2.pdf", dpi=600)
        plt.savefig("boxplotlim_2_corr.pdf", dpi=600)
        plt.savefig("boxplotlim_2_corr.eps", dpi=600)
        #plt.savefig("boxplotlim2.tex", dpi=600)
   
    plt.show()
    exit()

#########################################################
#########################################################

# # fft figure
# plt.figure()
# XS = np.fft.rfft(xs)
# freqs = np.fft.fftfreq(len(xs), dt)
# freqs[len(XS)-1] *= -1 # nyquist returned for negative freq
# freqs = freqs[0:len(XS)]
# 
# plt.plot(freqs[0:len(XS)], XS.real)
# plt.plot(freqs[0:len(XS)], XS.imag)
# 
# # truncate by T X
# plt.figure()
# dv = freqs[1] - freqs[0]
# iLw = vL / dv
# tr = int(len(freqs) / 35.)
# plt.plot(freqs[iLw-tr:iLw + tr], XS.real[iLw-tr:iLw+tr], '-x')
# plt.plot(freqs[iLw-tr:iLw + tr], XS.imag[iLw-tr:iLw+tr], '-x')
# plt.twinx()
# plt.plot(freqs[iLw-tr:iLw + tr], np.angle(XS)[iLw-tr:iLw+tr], '-x')
# 
# print len(XS[iLw-tr:iLw+tr])
# print len(t[iisamp::irsamp])
# 
# plt.show()