From d9beb94422ed5efce5901b5456cb10e1b13be1a6 Mon Sep 17 00:00:00 2001 From: Oleg Petrov Date: Wed, 30 Sep 2015 11:09:36 +0000 Subject: [PATCH] New experimetns added --- Scripts/EXSY/op_noesy_exp.py | 176 +++++++++++++++ Scripts/EXSY/op_noesy_res.py | 157 +++++++++++++ Scripts/EXSY_2H/op_exsy2h_exp.py | 194 ++++++++++++++++ Scripts/EXSY_2H/op_exsy2h_res.py | 170 ++++++++++++++ .../op_spinal4pulses_exp.py | 176 +++++++++++++++ .../op_spinal4pulses_res.py | 210 ++++++++++++++++++ .../op_sgse16_exp.py | 183 +++++++++++++++ .../op_sgse16_res.py | 201 +++++++++++++++++ .../op_zeeman4pulses_exp.py | 173 +++++++++++++++ .../op_zeeman4pulses_res.py | 210 ++++++++++++++++++ 10 files changed, 1850 insertions(+) create mode 100644 Scripts/EXSY/op_noesy_exp.py create mode 100644 Scripts/EXSY/op_noesy_res.py create mode 100644 Scripts/EXSY_2H/op_exsy2h_exp.py create mode 100644 Scripts/EXSY_2H/op_exsy2h_res.py create mode 100644 Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_exp.py create mode 100644 Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_res.py create mode 100644 Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_exp.py create mode 100644 Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_res.py create mode 100644 Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_exp.py create mode 100644 Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_res.py diff --git a/Scripts/EXSY/op_noesy_exp.py b/Scripts/EXSY/op_noesy_exp.py new file mode 100644 index 0000000..31e94a8 --- /dev/null +++ b/Scripts/EXSY/op_noesy_exp.py @@ -0,0 +1,176 @@ +# -*- coding: iso-8859-1 -*- + +TXEnableDelay = 2e-6 +TXEnableValue = 0b0001 # TTL line enabling RF amplifier (bit 0) +TXPulseValue = 0b0010 # TTL line triggering RF pulses (bit 1) +ADCSensitivity = 2 # voltage span for ADC + +def experiment(): # 2D NOESY experiment, using States-TPPI technique for quadrature detection in F1 + + # States-TPPI technique achieves two effects for an indirect dimension F1: + # (1) signal frequency discrimination and (2) displacement of the unmodulated + # artefact signal from an inconvenient location in the middle of spectrum to the edge. + # (1) is achieved by recording two data sets at each t1 point - with orthogonal phases + # of the preparation pulse and same receiver phase - and storing them in separate memory + # locations. These two fid measurements yield one complex data point in F1. + # (2) by inverting phase of the preparation pulse and the receiver each time when t1 is + # incremented (that is for subsequent complex points). Therefore, the artefact signal + # becomes modulated at the Nyquist frequency and appears in the spectrum at F1=ħSW/2 Hz + # instead of 0 Hz, where SW is spectral width. [http://nmrwiki.org] + + # set up acqusition parameters: + pars = {} + pars['P90'] = 1.65e-6 # 90-degree pulse length (s) + pars['SF'] = 338.7e6 # spectrometer frequency (Hz) + pars['O1'] = -57.0e3 # offset from SF (Hz) + pars['SW'] = 150e3 # spectral width (Hz) + pars['SI1'] = 32 # number of (complex) data points in F1 (2D) + pars['SI2'] = 128 # number of (complex) data points in F2 + pars['D8'] = 100e-6 # mixing time, tm (s) + pars['NS'] = 16 # number of scans + pars['DS'] = 0 # number of dummy scans + pars['RD'] = 2.5 # delay between scans (s) + pars['DEAD1'] = 4e-6 # receiver dead time (s) + pars['PHA'] = 150 # receiver reference phase (degree) + pars['DATADIR'] = '/home/fprak/' # data directory + pars['OUTFILE'] = 'test' # output file name + + # specify a variable parameter (optional): + pars['VAR_PAR'] = 'D8' # variable parameter name (a string) + start = 10.e-6 # starting value + stop = 1000e-6 # end value + steps = 3 # number of values + log_scale = True # log scale flag + stag_range = False # staggered range flag + + # check parameters for safety: + if pars['PHA'] < 0: + pars['PHA'] = 360 + pars['PHA'] + + if pars['P90'] > 20e-6: + raise Exception("Pulse too long!!!") + + var_key = pars.get('VAR_PAR') + if var_key == 'P90' and (start > 20e-6 or stop > 20e-6): + raise Exception("Pulse too long!!!") + + if pars['NS']%8 != 0: + pars['NS'] = (pars['NS'] / 8 + 1) * 8 + print 'Number of scans changed to', pars['NS'], 'due to phase cycling' + + # start the experiment: + if var_key: + # this is an arrayed experiment: + if log_scale: + array = log_range(start,stop,steps) + else: + array = lin_range(start,stop,steps) + + if stag_range: + array = staggered_range(array, size=2) + + # estimate the experiment time: + if var_key == 'D8': + seconds = (sum(array) + (.5*pars['SI1']/pars['SW'] + pars['RD']) * steps) * (pars['NS'] + pars['DS']) * 2*pars['SI1'] + elif var_key == 'RD': + seconds = (sum(array) + (.5*pars['SI1']/pars['SW'] + pars['D8']) * steps) * (pars['NS'] + pars['DS']) * 2*pars['SI1'] + else: + seconds = (.5*pars['SI1']/pars['SW'] + pars['D8'] + pars['RD']) * steps * (pars['NS']+ pars['DS']) * 2*pars['SI1'] + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + + # loop for a variable parameter: + for index, pars[var_key] in enumerate(array): + print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\ + ' out of '+str(array.size)+': value = '+str(pars[var_key]) + # loop for accumulation and sampling the indirect dimension F1: + for run in xrange((pars['NS']+pars['DS'])*2*pars['SI1']): + yield noesyst_experiment(pars, run) + synchronize() + + else: + # estimate the experiment time: + seconds = (.5*pars['SI1']/pars['SW'] + pars['D8'] + pars['RD']) * (pars['NS']+ pars['DS']) * 2*pars['SI1'] + print 'sec ', seconds + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for accumulation and sampling the indirect dimension F1: + for run in xrange((pars['NS']+pars['DS'])*2*pars['SI1']): + yield noesyst_experiment(pars, run) + + +# the pulse program: + +def noesyst_experiment(pars, run): + e=Experiment() + + dummy_scans = pars.get('DS') + if dummy_scans: + run -= dummy_scans + + # phase lists (M.H.Levitt 'Spin Dynamics', 2nd edition, p.530): + pars['PH1'] = [ 0, 180, 0, 180, 0, 180, 0, 180] + pars['PH3'] = [180, 180, 180, 180, 180, 180, 180, 180] + pars['PH4'] = [ 0, 0, 90, 90, 180, 180, 270, 270] + pars['PH2'] = [ 0, 180, 90, 270, 180, 0, 270, 90] # receiver + + # read in variables: + P90 = pars['P90'] + SF = pars['SF'] + O1 = pars['O1'] + RD = pars['RD'] + DEAD1 = pars['DEAD1'] + D8 = pars['D8'] + NS = pars['NS'] + PH1 = pars['PH1'][run%len(pars['PH1'])] + PH3 = pars['PH3'][run%len(pars['PH3'])] + PH4 = pars['PH4'][run%len(pars['PH4'])] + PH2 = pars['PH2'][run%len(pars['PH2'])] + PHA = pars['PHA'] + + # F1 sampling parameters: + IN0 = 1./pars['SW'] # t1 increment + + # the States-TPPI bit: + PH1-= (run/(1*NS))%4*90 # PH1 changes by 90-deg. after every 1*NS scans + D0 = (run/(2*NS)) *IN0 # t1 increases by IN0 after every 2*NS scans + + # F2 sampling parameters: + SI2 = pars['SI2'] + SW2 = pars['SW'] + while SW2 <= 10e6 and SI2 < 256*1024: + SI2 *= 2 + SW2 *= 2 + + # run the pulse sequence: + e.wait(RD) # delay between scans + e.set_frequency(SF+O1, phase=PH1) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(D0) # incremented delay t1 + e.set_phase(PH3) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(D8) # mixing time + e.set_phase(PH4) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + e.set_phase(PHA) # set phase for receiver + e.wait(DEAD1) # wait for coil ringdown + e.record(SI2, SW2, sensitivity=ADCSensitivity) # acquire signal + + # write experiment attributes: + for key in pars.keys(): + e.set_description(key, pars[key]) # acqusition parameters + e.set_description('run', run) # current scan + e.set_description('rec_phase', -PH2) # current receiver phase + + return e \ No newline at end of file diff --git a/Scripts/EXSY/op_noesy_res.py b/Scripts/EXSY/op_noesy_res.py new file mode 100644 index 0000000..f39c79a --- /dev/null +++ b/Scripts/EXSY/op_noesy_res.py @@ -0,0 +1,157 @@ +# -*- coding: iso-8859-1 -*- + +from numpy import * +from scipy.signal import * +from scipy.optimize import * +from os import path, rename + +def result(): + + measurement = MeasurementResult('Magnetization') + + measurement_range = [0.0, 10e-6] + measurement_ranging = False + + suffix = '' # output file name's suffix and... + counter = 0 # counter for arrayed 2D experiments + + # loop over the incoming results: + for timesignal in results: + if not isinstance(timesignal,ADC_Result): + continue + + # read experiment parameters: + pars = timesignal.get_description_dictionary() + + # ---------------- digital filter ------------------ + + # get actual sampling rate of timesignal: + sampling_rate = timesignal.get_sampling_rate() + + # get user-defined spectrum width: + spec_width = pars['SW'] + + # specify cutoff frequency, in relative units: + cutoff = spec_width / sampling_rate + + if cutoff < 1: # otherwise no filter applied + + # number of filter's coefficients: + numtaps = 29 + + # use firwin to create a lowpass FIR filter: + fir_coeff = firwin(numtaps, cutoff) + + # downsize x according to user-defined spectral window: + skip = int(sampling_rate / spec_width) + timesignal.x = timesignal.x[::skip] + + for i in range(2): + # apply the filter to ith channel: + timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i]) + + # zeroize first N-1 "corrupted" samples: + timesignal.y[i][:numtaps-1] = 0.0 + + # circular left shift of y: + timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1)) + + # downsize y to user-defined number of samples (SI): + timesignal.y[i] = timesignal.y[i][::skip] + + # update the sampling_rate attribute of the signal's: + timesignal.set_sampling_rate(spec_width) + + # ---------------------------------------------------- + + # rotate timesignal according to current receiver's phase: + timesignal.phase(pars['rec_phase']) + + # provide timesignal to the display tab: + data['Current scan'] = timesignal + + # accumulate... + if not locals().get('accu'): + accu = Accumulation() + + # skip dummy scans, if any: + if pars['run'] < 0: continue + + # add up: + accu += timesignal + + # provide accumulation to the display tab: + data['Accumulation'] = accu + + # check whether accumulation is complete: + # ---------------------------------------------------------------------------------------- + # Note that State-TPPI technique implies recording two data sets at each t1 point. + # Cosine-modulated data are stored in record 1 as Re, and sine-modulated data are stored + # in record 2 as Im, totally 2*SI1 records. Henceforth, accu represents one such a record. + # ----------------------------------------------------------------------------------------- + if accu.n == pars['NS']: + + # make a copy: + fid = accu + 0 + + # compute the initial phase of FID: + phi0 = arctan2(fid.y[1][0], fid.y[0][0]) * 180 / pi + if not 'ref' in locals(): ref = phi0 + # print 'phi0 = ', phi0 + + # rotate FID to maximize y[0][0]: + #echo.phase(-phi0) + + # check whether it is an arrayed experiment: + var_key = pars.get('VAR_PAR') + if var_key: + # get variable parameter's value: + var_value = pars.get(var_key) + + # provide signal recorded with this var_value to the display tab: + data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu + + # measure signal intensity vs. var_value: + if measurement_ranging == True: + [start, stop] = accu.get_sampling_rate() * array(measurement_range) + measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) + + else: + measurement[var_value] = sum(accu.y[0][0:31]) + + # provide measurement to the display tab: + data[measurement.get_title()] = measurement + + # update the file name suffix: + counter2D = counter/(2*pars['SI1'])+1 + suffix = '_' + str(counter2D) + counter += 1 + + # save accu if required: + outfile = pars.get('OUTFILE') + if outfile: + datadir = pars.get('DATADIR') + + # write raw data in Tecmag format: + filename = datadir+outfile+suffix+'.tnt' + accu.write_to_tecmag(filename,\ + nrecords=2*pars['SI1'],\ + frequency=pars['SW']+pars['O1']) + + # write parameters in a text file: + filename = datadir+outfile+suffix+'.par' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.par') + + fileobject = open(filename, 'w') + for key in sorted(pars.iterkeys()): + if key=='run': continue + if key=='rec_phase': continue + fileobject.write('%s%s%s'%(key,'=', pars[key])) + fileobject.write('\n') + fileobject.close() + + # reset accumulation: + del accu + +pass \ No newline at end of file diff --git a/Scripts/EXSY_2H/op_exsy2h_exp.py b/Scripts/EXSY_2H/op_exsy2h_exp.py new file mode 100644 index 0000000..7f57d08 --- /dev/null +++ b/Scripts/EXSY_2H/op_exsy2h_exp.py @@ -0,0 +1,194 @@ +# -*- coding: iso-8859-1 -*- + +TXEnableDelay = 0.5e-6 +TXEnableValue = 0b0001 # TTL line blanking RF amplifier (bit 0) +TXPulseValue = 0b0010 # TTL line triggering RF pulses (bit 1) +ADCSensitivity = 1 # voltage span for ADC + +def experiment(): # 2H Exchange Spectroscopy (2H EXSY) experiment [JMR 79, 269-290 (1988)] + + # Cosine and sine modulated signals are acquired sequentially by switching + # between Zeeman order and spin-alignment phase lists. The signals are + # processed into a pure absorption mode 2D spectrum according to scheme by + # [Bluemich, Schmidt, and Spiess, JMR 79, 269-290 (1988)]. Prior to writing + # in a file (Tecmag), the sine-modulated signal is rotated by 90°, thus + # enabling 2D processing via a regular States algorithm with NMRnotebook + # or like NMR software. + + # set up acquisition parameters: + pars = {} + pars['P90'] = 2.7e-6 # 90°-pulse length (s) + pars['SF'] = 46.140e6 # spectrometer frequency (Hz) + pars['O1'] = 1000 # offset from SF (Hz) + pars['SW'] = 125e3 # spectral window (Hz) + pars['SI1'] = 80 # number of (complex) data points in F1 (2nd dimension) + pars['SI2'] = 1*256 # number of (complex) data points in F2 + pars['D3'] = 10e-6 # position of refocusing 90°-pulse, Delta (s) + pars['D4'] = 2e-6 # pre-aquisition delay (s) + pars['D8'] = 3e-3 # mixing time, tm (s) + pars['NS'] = 512 # number of scans + pars['DS'] = 0 # number of dummy scans + pars['RD'] = 0.2 # delay between scans (s) + pars['PHA'] = 65 # receiver reference phase (degree) + pars['DATADIR'] = '/home/mathilda/Desktop/Oleg/temp/' # data directory + pars['OUTFILE'] = 'dso_320K' # output file name + + # specify a variable parameter (optional): + pars['VAR_PAR'] = None # variable parameter's name (a string) + start = 80 # starting value + stop = 128 # end value + steps = 2 # number of values + log_scale = False # log scale flag + stag_range = False # staggered range flag + + # check parameters for safety: + if pars['PHA'] < 0: + pars['PHA'] = 360 + pars['PHA'] + + if pars['P90'] > 20e-6: + raise Exception("Pulse too long!!!") + + # check whether a variable parameter is named: + var_key = pars.get('VAR_PAR') + if var_key == 'P90' and (start > 20e-6 or stop > 20e-6): + raise Exception("Pulse too long!!!") + + if pars['NS']%16 != 0: + pars['NS'] = int(round(pars['NS'] / 32) + 1) * 32 + print 'Number of scans changed to ', pars['NS'], ' due to phase cycling' + + # start the experiment: + if var_key: + # this is an arrayed experiment: + if log_scale: + array = log_range(start,stop,steps) + else: + array = lin_range(start,stop,steps) + + if stag_range: + array = staggered_range(array, size = 2) + + # estimate the experiment time: + if var_key == 'D8': + seconds = (sum(array) + (.5*pars['SI1']/pars['SW'] + pars['RD']) * steps) * (pars['NS'] + pars['DS']) * 2*pars['SI1'] + elif var_key == 'RD': + seconds = (sum(array) + (.5*pars['SI1']/pars['SW'] + pars['D8']) * steps) * (pars['NS'] + pars['DS']) * 2*pars['SI1'] + else: + seconds = (.5*pars['SI1']/pars['SW'] + pars['D8'] + pars['RD']) * steps * (pars['NS']+ pars['DS']) * 2*pars['SI1'] + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for a variable parameter: + for index, pars[var_key] in enumerate(array): + print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\ + ' out of '+str(array.size)+': value = '+str(pars[var_key]) + # loop for accumulation and sampling the indirect dimension F1: + for run in xrange((pars['NS']+pars['DS'])*2*pars['SI1']): + yield exsy2h_experiment(pars, run) + synchronize() + + else: + # estimate the experiment time: + seconds = (.5*pars['SI1']/pars['SW'] + pars['D8'] + pars['RD']) * (pars['NS']+ pars['DS']) * 2*pars['SI1'] + print 'sec ', seconds + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for accumulation and sampling the indirect dimension F1: + for run in xrange((pars['NS']+pars['DS'])*2*pars['SI1']): + yield exsy2h_experiment(pars, run) + + +# the pulse program: + +def exsy2h_experiment(pars, run): + e=Experiment() + + dummy_scans = pars.get('DS') + if dummy_scans: + run -= dummy_scans + + pars['PROG'] = 'exsy2h_experiment' + + # 8-step phase cycle (1-21) modifided to deal with T1-recovery and extended for Re/Im imbalance) + pars['PH1'] = [0, 270, 0, 270, 180, 90, 180, 90] # 1st pulse (90°) + pars['PH3'] = [0, 90, 0, 90, 0, 90, 0, 90] # 2nd pulse (57.4°) + pars['PH4'] = [0, 0, 180, 180, 270, 270, 90, 90] # 3rd pulse (57.4°) + pars['PH5'] = [90, 90, 90, 90, 180, 180, 180, 180] # 4th pulse (90°) + pars['PH2'] = [0, 180, 180, 0, 90, 270, 270, 90] # receiver + + # read in variables: + P90 = pars['P90'] + P1 = pars['P90']*(54.7/90) + SF = pars['SF'] + O1 = pars['O1'] + RD = pars['RD'] + D4 = pars['D4'] + D8 = pars['D8'] + D3 = pars['D3'] + NS = pars['NS'] + PH1 = pars['PH1'][run%len(pars['PH1'])] + PH3 = pars['PH3'][run%len(pars['PH3'])] + PH4 = pars['PH4'][run%len(pars['PH4'])] + PH5 = pars['PH5'][run%len(pars['PH5'])] + PH2 = pars['PH2'][run%len(pars['PH2'])] + PHA = pars['PHA'] + + # this is a part of phase cycling: + PH5 += (run/len(pars['PH5']))%2*180 + PH1 += (run/len(pars['PH5']))%2*180 + PH2 += (run/len(pars['PH5']))%2*180 + + # F1 sampling parameters: + IN0 = 1./pars['SW'] # t1 increment + + # F1 sampling scheme: + PH3+= (run/(1*NS))%4*90 # phases are upgraded after every NS scans + PH4+= (run/(1*NS))%4*90 + D0 = (run/(2*NS)) *IN0 # t1 is incremented after every 2*NS scans + + # F2 sampling parameters: + SI2 = pars['SI2'] + SW2 = pars['SW'] + while SW2 <= 10e6 and SI2 < 256*1024: + SI2 *= 2 + SW2 *= 2 + + # run the pulse sequence: + e.wait(RD) # relaxation delay between scans + e.set_frequency(SF+O1, phase=PH1) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90°-pulse + + e.wait(D0) # incremented delay, t1 + e.set_phase(PH3) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P1, value=TXEnableValue|TXPulseValue) # 54.7°-pulse + + e.wait(D8) # mixing time, tm + e.set_phase(PH4) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P1, value=TXEnableValue|TXPulseValue) # 54.7°-pulse + + e.wait(D3) # refocusing delay, Delta + + e.set_phase(PH5) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90°-pulse + + e.wait(TXEnableDelay) + e.set_phase(PHA) + e.wait(D3+D4) # pre-aquisition delay + e.record(SI2, SW2, sensitivity=ADCSensitivity) # acquisition + + # write experiment parameters: + for key in pars.keys(): + e.set_description(key, pars[key]) # acquisition parameters + e.set_description('run', run) # current scan + e.set_description('rec_phase', -PH2) # current receiver phase + + return e \ No newline at end of file diff --git a/Scripts/EXSY_2H/op_exsy2h_res.py b/Scripts/EXSY_2H/op_exsy2h_res.py new file mode 100644 index 0000000..8afe909 --- /dev/null +++ b/Scripts/EXSY_2H/op_exsy2h_res.py @@ -0,0 +1,170 @@ +# -*- coding: iso-8859-1 -*- + +from numpy import * +from scipy.signal import * +from scipy.optimize import * +from scipy.fftpack import fft, ifft +from os import path, rename + +def result(): + + measurement = MeasurementResult('Magnetization') + + measurement_range = [0.0, 10e-6] + measurement_ranging = False + + suffix = '' # output file name's suffix and... + counter = 0 # counter for arrayed 2D experiments + # npts = 0 + + # loop over the incoming results: + for timesignal in results: + if not isinstance(timesignal,ADC_Result): + continue + + # read experiment parameters: + pars = timesignal.get_description_dictionary() + + # ---------------- digital filter ------------------ + + # get actual sampling rate of timesignal: + sampling_rate = timesignal.get_sampling_rate() + + # get user-defined spectrum width: + spec_width = pars['SW'] + + # specify cutoff frequency, in relative units: + cutoff = spec_width / sampling_rate + + if cutoff < 1: # otherwise no filter applied + + # number of filter's coefficients: + numtaps = 29 + + # use firwin to create a lowpass FIR filter: + fir_coeff = firwin(numtaps, cutoff) + + # downsize x according to user-defined spectral window: + skip = int(sampling_rate / spec_width) + timesignal.x = timesignal.x[::skip] + + for i in range(2): + # apply the filter to ith channel: + timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i]) + + # zeroize first N-1 "corrupted" samples: + timesignal.y[i][:numtaps-1] = 0.0 + + # circular left shift of y: + timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1)) + + # downsize y to user-defined number of samples (SI): + timesignal.y[i] = timesignal.y[i][::skip] + + # update the sampling_rate attribute of the signal's: + timesignal.set_sampling_rate(spec_width) + + # ---------------------------------------------------- + + # rotate timesignal according to current receiver's phase: + timesignal.phase(pars['rec_phase']) + + # provide timesignal to the display tab: + data['Current scan'] = timesignal + + # accumulate... + if not locals().get('accu'): + accu = Accumulation() + + # skip dummy scans, if any: + if pars['run'] < 0: continue + + # add up: + accu += timesignal + + # provide accumulation to the display tab: + data['Accumulation'] = accu + + # check whether accumulation is complete: + # ------------------------------------------------------------------------------------ + # The hypercomplex technique implies recording two data sets for each t1 value. + # One dataset (a cosine-modulated signal) is stored in odd records as Re, while + # the other dataset (a sine-modulated signal) in even records as Im, totally 2*SI1 + # records. Henceforth, accu represents one such a record. + # ------------------------------------------------------------------------------------ + if accu.n == pars['NS']: + # compute the initial phase of FID: + phi0 = arctan2(accu.y[1][0], accu.y[0][0]) * 180 / pi + if not 'ref' in locals(): ref = phi0 + print 'phi0 = ', phi0 + + # rotate every other record by 90° so that States algorithm is applicable: + rec = (accu.job_id/accu.n)%(2*pars['SI1']) + 1 + if rec%2 == 0: + accu.phase(90) + coeff = 1.5 + accu.y[0] *= coeff # XY-balancing + accu.y[1] *= coeff + + else: # baseline correction ))))) + tmp = fft(accu.y[0]+1j*accu.y[1]) + [start, stop] = len(accu.y[0])*array([0.4, 0.6]) + tmp -= mean(tmp.real[start:stop]) + tmp = ifft(tmp) + accu.y[0] = tmp.real + accu.y[1] = tmp.imag + + # check whether it is an arrayed experiment: + var_key = pars.get('VAR_PAR') + if var_key: + # get variable parameter's value: + var_value = pars.get(var_key) + + # provide signal recorded with this var_value to the display tab: + data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu + + # measure signal intensity vs. var_value: + if measurement_ranging == True: + [start, stop] = accu.get_sampling_rate() * array(measurement_range) + measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) + + else: + # measurement[var_value] = sum(accu.y[0][0:31]) + measurement[var_value+counter*1e-6] = sum(accu.y[0][0:31]) + + # provide measurement to the display tab: + data[measurement.get_title()] = measurement + + # update the file name suffix: + counter2D = counter/(2*pars['SI1'])+1 + suffix = '_' + str(counter2D) + counter += 1 + + # save accu if required: + outfile = pars.get('OUTFILE') + if outfile: + datadir = pars.get('DATADIR') + + # write raw data in Tecmag format: + filename = datadir+outfile+suffix+'.tnt' + accu.write_to_tecmag(filename,\ + nrecords=2*pars['SI1'],\ + frequency=pars['SW']+pars['O1']) + + # write parameters in a text file: + filename = datadir+outfile+suffix+'.par' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.par') + + fileobject = open(filename, 'w') + for key in sorted(pars.iterkeys()): + if key=='run': continue + if key=='rec_phase': continue + fileobject.write('%s%s%s'%(key,'=', pars[key])) + fileobject.write('\n') + fileobject.close() + + # reset accumulation: + del accu + +pass \ No newline at end of file diff --git a/Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_exp.py b/Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_exp.py new file mode 100644 index 0000000..0fc8dbf --- /dev/null +++ b/Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_exp.py @@ -0,0 +1,176 @@ +# -*- coding: iso-8859-1 -*- + +TXEnableDelay = 1e-6 +TXEnableValue = 0b0001 # TTL line blanking RF amplifier (bit 0) +TXPulseValue = 0b0010 # TTL line triggering RF pulses (bit 1) +ADCSensitivity = 1 # voltage span for ADC + +def experiment(): # Jeener-Broekaert echo sequence (a.k.a. spin-alignment) + + # set up acquisition parameters: + pars = {} + pars['P90'] = 4.2e-6 # 90-degree pulse length (s) + pars['SF'] = 46.7e6 # spectrometer frequency (Hz) + pars['O1'] = -30e3 # offset from SF (Hz) + pars['SW'] = 10e6 # spectral window (Hz) + pars['SI'] = 1*1024 # number of acquisition points + pars['NS'] = 16 # number of scans + pars['DS'] = 0 # number of dummy scans + pars['RD'] = 0.5 # delay between scans (s) + pars['D1'] = 30e-6 # delay after first pulse, or tp (s) + pars['D2'] = 100e-6 # delay after second pulse, or tm (s) + pars['D3'] = 20e-6 # refocusing pulse delay (s) + pars['PHA'] = 30 # receiver reference phase (degree) + pars['DATADIR'] = '/home/fprak/Students/' # data directory + pars['OUTFILE'] = None # output file name + + # specify a variable parameter (optional): + pars['VAR_PAR'] = 'D2' # variable parameter name (a string) + start = 30e-6 # starting value + stop = 2e-0 # end value + steps = 24 # number of values + log_scale = True # log scale flag + stag_range = False # staggered range flag + + # check parameters for safety: + if pars['PHA'] < 0: + pars['PHA'] = 360 + pars['PHA'] + + if pars['P90'] > 20e-6: + raise Exception("Pulse too long!!!") + + # check whether a variable parameter is named: + var_key = pars.get('VAR_PAR') + if var_key == 'P90' and (start > 20e-6 or stop > 20e-6): + raise Exception("Pulse too long!!!") + + if pars['NS']%16 != 0: + pars['NS'] = int(round(pars['NS'] / 16) + 1) * 16 + print 'Number of scans changed to ', pars['NS'], ' due to phase cycling' + + # start the experiment: + if var_key: + # this is an arrayed experiment: + if log_scale: + array = log_range(start,stop,steps) + else: + array = lin_range(start,stop,steps) + + if stag_range: + array = staggered_range(array, size = 2) + + # estimate the experiment time: + if var_key == 'D1': + seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS']) + elif var_key == 'D2': + seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS']) + elif var_key == 'RD': + seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS']) + else: + seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS']) + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for a variable parameter: + for index, pars[var_key] in enumerate(array): + print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\ + ' out of '+str(array.size)+': value = '+str(pars[var_key]) + # loop for accumulation: + for run in xrange(pars['NS']+pars['DS']): + yield spinal4pulses_experiment(pars, run) + synchronize() + + else: + # estimate the experiment time: + seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS']) + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for accumulation: + for run in xrange(pars['NS']+pars['DS']): + yield spinal4pulses_experiment(pars, run) + + +# the pulse program: + +def spinal4pulses_experiment(pars, run): + e=Experiment() + + dummy_scans = pars.get('DS') + if dummy_scans: + run -= dummy_scans + + pars['PROG'] = 'spinal4pulses_experiment' + + # 8-step phase cycle (1-14 modifided to deal with T1-recovery and extended for Re/Im imbalance) + pars['PH1'] = [0, 270, 0, 270, 180, 90, 180, 90] # 1st (90-degree) pulse + pars['PH3'] = [90, 180, 90, 180, 90, 180, 90, 180] # 2nd (90-degree) pulse + pars['PH4'] = [90, 90, 270, 270, 0, 0, 180, 180] # 3rd (90-degree) pulse + pars['PH5'] = [90, 90, 90, 90, 180, 180, 180, 180] # 3rd (90-degree) pulse + pars['PH2'] = [0, 180, 180, 0, 90, 270, 270, 90] # receiver + + # read in variables: + P90 = pars['P90'] + P45 = pars['P90']*0.5 + P1 = pars['P90']*0.5 + SF = pars['SF'] + O1 = pars['O1'] + RD = pars['RD'] + D1 = pars['D1'] + D2 = pars['D2'] + D3 = pars['D3'] + PH1 = pars['PH1'][run%len(pars['PH1'])] + PH3 = pars['PH3'][run%len(pars['PH3'])] + PH4 = pars['PH4'][run%len(pars['PH4'])] + PH5 = pars['PH5'][run%len(pars['PH5'])] + PH2 = pars['PH2'][run%len(pars['PH2'])] + PHA = pars['PHA'] + + if (run/8)%2 != 0: + PH5 += 180 + + # set sampling parameters: + SI = pars['SI'] + SW = pars['SW'] + while SW <= 10e6 and SI < 256*1024: + SI *= 2 + SW *= 2 + + # run the pulse sequence: + e.wait(RD) # relaxation delay between scans + e.set_frequency(SF+O1, phase=PH1) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(D1-P90/2-P45/2-TXEnableDelay) # 'short tau' + e.set_phase(PH3) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P45, value=TXEnableValue|TXPulseValue) # 45-degree pulse + + e.wait(D2-P45/2-P1/2-TXEnableDelay) # 'long tau' + e.set_phase(PH4) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P1, value=TXEnableValue|TXPulseValue) # 45-degree pulse + + e.wait(D3-P1/2-P90/2-TXEnableDelay) # 'delta' + + e.set_phase(PH5) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(TXEnableDelay) + e.set_phase(PHA) + e.wait(D3+D1-P90/2-TXEnableDelay) # wait for echo + e.record(SI, SW, sensitivity=ADCSensitivity) # acquisition + + # write experiment parameters: + for key in pars.keys(): + e.set_description(key, pars[key]) # acquisition parameters + e.set_description('run', run) # current scan + e.set_description('rec_phase', -PH2) # current receiver phase + + return e \ No newline at end of file diff --git a/Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_res.py b/Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_res.py new file mode 100644 index 0000000..5dc8a3c --- /dev/null +++ b/Scripts/Spin_Alignment_Four_Pulses/op_spinal4pulses_res.py @@ -0,0 +1,210 @@ +# -*- coding: iso-8859-1 -*- + +from numpy import * +from scipy.signal import * +from scipy.optimize import * +from os import path, rename + +def result(): + + measurement = MeasurementResult('Magnetization') + + measurement_range = [0.0, 10e-6] + measurement_ranging = False + + suffix = '' # output file name's suffix and... + counter = 1 # counter for arrayed experiments + var_key = '' # variable parameter name + + # loop over the incoming results: + for timesignal in results: + if not isinstance(timesignal,ADC_Result): + continue + + # read experiment parameters: + pars = timesignal.get_description_dictionary() + + # ---------------- digital filter ------------------ + + # get actual sampling rate of timesignal: + sampling_rate = timesignal.get_sampling_rate() + + # get user-defined spectrum width: + spec_width = pars['SW'] + + # specify cutoff frequency, in relative units: + cutoff = spec_width / sampling_rate + + if cutoff < 1: # no filter applied otherwise + + # number of filter's coefficients: + numtaps = 29 + + # use firwin to create a lowpass FIR filter: + fir_coeff = firwin(numtaps, cutoff) + + # downsize x according to user-defined spectral window: + skip = int(sampling_rate / spec_width) + timesignal.x = timesignal.x[::skip] + + for i in range(2): + # apply the filter to ith channel: + timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i]) + + # zeroize first N-1 "corrupted" samples: + timesignal.y[i][:numtaps-1] = 0.0 + + # circular left shift of y: + timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1)) + + # downsize y to user-defined number of samples (SI): + timesignal.y[i] = timesignal.y[i][::skip] + + # update the sampling_rate attribute of the signal's: + timesignal.set_sampling_rate(spec_width) + + # ---------------------------------------------------- + + # rotate timesignal according to current receiver's phase: + timesignal.phase(pars['rec_phase']) + + # provide timesignal to the display tab: + data['Current scan'] = timesignal + + # accumulate... + if not locals().get('accu'): + accu = Accumulation() + + # skip dummy scans, if any: + if pars['run'] < 0: continue + + # add up: + accu += timesignal + + # provide accumulation to the display tab: + data['Accumulation'] = accu + + # check how many scans are done: + if accu.n == pars['NS']: # accumulation is complete + + # make a copy: + echo = accu + 0 + + # compute the initial phase of FID: + phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi + if not 'ref' in locals(): ref = phi0 + print 'phi0 = ', phi0 + + # rotate FID to maximize y[0][0]: + #echo.phase(-phi0) + + # do FFT: + echo.exp_window(line_broadening=10) + spectrum = echo.fft(samples=2*pars['SI']) + + # try zero-order phase correction: + spectrum.phase(-phi0) + + # provide spectrum to the display tab: + data['Spectrum'] = spectrum + + # check whether it is an arrayed experiment: + var_key = pars.get('VAR_PAR') + if var_key: + # get variable parameter's value: + var_value = pars.get(var_key) + + # provide signal recorded with this var_value to the display tab: + data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu + + # measure signal intensity vs. var_value: + if measurement_ranging == True: + [start, stop] = accu.get_sampling_rate() * array(measurement_range) + measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) + + else: + measurement[var_value] = sum(accu.y[0][0:31]) + + # provide measurement to the display tab: + data[measurement.get_title()] = measurement + + # update the file name suffix: + suffix = '_' + str(counter) + counter += 1 + + # save accu if required: + outfile = pars.get('OUTFILE') + if outfile: + datadir = pars.get('DATADIR') + + # write raw data in Simpson format: + filename = datadir+outfile+suffix+'.dat' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.dat') + accu.write_to_simpson(filename) + + # write raw data in Tecmag format: + # filename = datadir+outfile+'.tnt' + # accu.write_to_tecmag(filename, nrecords=20) + + # write parameters in a text file: + filename = datadir+outfile+suffix+'.par' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.par') + + fileobject = open(filename, 'w') + for key in sorted(pars.iterkeys()): + if key=='run': continue + if key=='rec_phase': continue + fileobject.write('%s%s%s'%(key,'=', pars[key])) + fileobject.write('\n') + fileobject.close() + + # reset accumulation: + del accu + + + if var_key == 'D2': + # KWW fit: + xdata = measurement.get_xdata() + ydata = measurement.get_ydata() + [amplitude, rate, beta] = fitfunc(xdata, ydata) + print '%s%02g' % ('Amplitude = ', amplitude) + print '%s%02g' % ('T2 [s] = ', 1./rate) + print '%s%02g' % ('Beta = ', beta) + + # update display for the fit: + measurement.y = func([amplitude, rate, beta], xdata) + data[measurement.get_title()] = measurement + +# the fitting procedure: +def fitfunc(xdata, ydata): + + # initialize variable parameters: + try: + # solve Az = b: + A = array((ones(xdata.size/2), xdata[0:xdata.size/2])) + b = log(abs(ydata[0:xdata.size/2])) + z = linalg.lstsq(transpose(A), b) + amplitude = exp(z[0][0]) + rate = -z[0][1] + except: + amplitude = abs(ydata[0]) + rate = 1./(xdata[-1] - xdata[0]) + beta = 1 + p0 = [amplitude, rate, beta] + + # run least-squares optimization: + plsq = leastsq(residuals, p0, args=(xdata, ydata)) + + return plsq[0] # best-fit parameters + +def residuals(p, xdata, ydata): + return ydata - func(p, xdata) + +# here is the function to fit: +def func(p, xdata): + return p[0]*exp(-(p[1]*xdata)**p[2]) + + +pass \ No newline at end of file diff --git a/Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_exp.py b/Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_exp.py new file mode 100644 index 0000000..33f32cd --- /dev/null +++ b/Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_exp.py @@ -0,0 +1,183 @@ +# -*- coding: iso-8859-1 -*- + +TXEnableDelay = 2e-6 +TXEnableValue = 0b0001 # TTL line blanking RF amplifier (bit 0) +TXPulseValue = 0b0010 # TTL line triggering RF pulses (bit 1) +ADCSensitivity = 5 # voltage span for ADC + +def experiment(): # a diffusion editing sequence with stimulated echo and XY-16 detection + + # set up acqusition parameters: + pars = {} + pars['P90'] = 0.8e-6 # 90-degree pulse length (s) + pars['SF'] = 161.85e6 # spectrometer frequency (Hz) + pars['O1'] = 0e3 # offset from SF (Hz) + pars['NS'] = 16 # number of scans + pars['DS'] = 0 # number of dummy scans + pars['RD'] = 6 # delay between scans (s) + pars['D1'] = 20e-6 # delay after first STE pulse (s) + pars['D2'] = 30e-6 # delay after second STE pulse (s) + pars['D4'] = 2.5e-6 # pre-acquisition offset + pars['NECH'] = 128 # number of 180-degree pulses + pars['TAU'] = 50e-6 # half pulse period (s) + pars['PHA'] = -95 # receiver phase (degree) + pars['DATADIR'] = '/home/fprak/Desktop/test/' # data directory + pars['OUTFILE'] = None # output file name + + # specify a variable parameter (optional): + pars['VAR_PAR'] = 'D2' # variable parameter name (a string) + start = 20e-6 # starting value + stop = 3e-1 # end value + steps = 24 # number of values + log_scale = True # log scale flag + stag_range = False # staggered range flag + + # check parameters for safety: + if pars['PHA'] < 0: + pars['PHA'] = 360 + pars['PHA'] + + if pars['P90'] > 20e-6: + raise Exception("Pulse too long!!!") + + # check whether a variable parameter is named: + var_key = pars.get('VAR_PAR') + if var_key == 'P90' and (start > 20e-6 or stop > 20e-6): + raise Exception("Pulse too long!!!") + + if pars['NS']%16 != 0: + pars['NS'] = int(round(pars['NS'] / 16) + 1) * 16 + print 'Number of scans changed to ',pars['NS'],' due to phase cycling' + + # start the experiment: + if var_key: + # this is an arrayed experiment: + if log_scale: + array = log_range(start,stop,steps) + else: + array = lin_range(start,stop,steps) + + if stag_range: + array = staggered_range(array, size = 2) + + # estimate the experiment time: + if var_key == 'D1': + seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS']) + elif var_key == 'D2': + seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS']) + elif var_key == 'RD': + seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS']) + else: + seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS']) + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for a variable parameter: + for index, pars[var_key] in enumerate(array): + print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\ + ' out of '+str(array.size)+': value = '+str(pars[var_key]) + # loop for accumulation: + for run in xrange(pars['NS']+pars['DS']): + yield ste16_experiment(pars, run) + synchronize() + else: + # estimate the experiment time: + seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS']) + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for accumulation: + for run in xrange(pars['NS']+pars['DS']): + yield ste16_experiment(pars, run) + + +# the pulse program: + +def ste16_experiment(pars, run): + e=Experiment() + + dummy_scans = pars.get('DS') + if dummy_scans: + run -= dummy_scans + + pars['PROG'] = 'ste16_experiment' + + # phase lists [from J. Magn. Reson. 157, 31 (2002)]: + pars['PH1'] = [0, 180, 0, 180, 0, 180, 0, 180, 90, 270, 90, 270, 90, 270, 90, 270] # 1st 90-degree pulse + pars['PH3'] = [0, 0, 180, 180, 0, 0, 180, 180, 0, 0, 180, 180, 0, 0, 180, 180] # 2nd 90-degree pulse + pars['PH4'] = [0, 0, 0, 0, 180, 180, 180, 180, 0, 0, 0, 0, 180, 180, 180, 180] # 3nd 90-degree pulse + pars['PH5'] = [0, 0, 0, 0, 0, 0, 0, 0, 270, 270, 270, 270, 270, 270, 270, 270] # 180-degree pulses + pars['PH2'] = [0, 180, 180, 0, 180, 0, 0, 180, 270, 90, 90, 270, 90, 270, 270, 90] # receiver + + # XY-16 phase cycle + pars['PH5_XY16'] = [0, 90, 0, 90, 90, 0, 90, 0, 180, 270, 180, 270, 270, 180, 270, 180] # XY-16: 180-degree pulses + pars['PH2_XY16'] = [90, 270, 270, 90, 270, 270, 90, 90, 90, 270, 270, 90, 270, 270, 90, 90] # XY_16: receiver + + # read in variables: + P90 = pars['P90'] + P180 = pars['P90']*2 + SF = pars['SF'] + O1 = pars['O1'] + RD = pars['RD'] + D1 = pars['D1'] + D2 = pars['D2'] + D4 = pars['D4'] + TAU = pars['TAU'] + NECH = pars['NECH'] + PH1 = pars['PH1'][run%len(pars['PH1'])] + PH3 = pars['PH3'][run%len(pars['PH3'])] + PH4 = pars['PH4'][run%len(pars['PH4'])] + PH5 = pars['PH5'][run%len(pars['PH5'])] + PH2 = pars['PH2'][run%len(pars['PH2'])] + PH5_XY16 = [x+PH5 for x in pars['PH5_XY16']*(NECH/16)] + PH2_XY16 = [x+PH2 for x in pars['PH2_XY16']*(NECH/16)] + PHA = pars['PHA'] + + # set sampling parameters: + SI = 128 # number of echo samples + SW = 20e6 # sample rate + AQ = (SI+6)/SW # acquisition window + + if TAU < (P90+P180)/2+TXEnableDelay or TAU < (P180+TXEnableDelay+AQ)/2: + raise Exception('pulse period is too short!') + + if 2*TAU < P180+TXEnableDelay+SI/SW: + raise Exception('pulse period too short!') + + # run the pulse sequence: + e.wait(RD) # delay between scans + e.set_frequency(SF+O1, phase=PH1) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 1st 90-degree pulse + + e.wait(D1-P90/2-TXEnableDelay) # short delay + e.set_phase(PH3) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 2nd 90-degree pulse + + e.wait(D2-P90/2-TXEnableDelay) # long delay + e.set_phase(PH4) + + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 3rd 90-degree pulse + + e.wait(D1+TAU-P90/2-TXEnableDelay) # wait for first echo and tau + + for i in range(NECH): + e.set_phase(PH5_XY16[i]) # set phase for 180-degree pulse + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) # enable RF amplifier + e.ttl_pulse(P180, value=TXEnableValue|TXPulseValue) # apply 180-degree pulse + e.set_phase(PHA) # set phase for receiver + e.wait(TAU-(P180+TXEnableDelay+AQ)/2+D4) # pre-acquisition delay + e.record(SI, SW, timelength=AQ, sensitivity=ADCSensitivity) # echo acquisition + e.wait(TAU-(P180+TXEnableDelay+AQ)/2-D4) # post-acquisition delay + + # write experiment attributes: + for key in pars.keys(): + e.set_description(key, pars[key]) # acqusition parameters + e.set_description('run', run) # current scan + e.set_description('rec_phase', [-1*x for x in PH2_XY16]) # current phase list for data routing + + return e \ No newline at end of file diff --git a/Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_res.py b/Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_res.py new file mode 100644 index 0000000..cd622b6 --- /dev/null +++ b/Scripts/Steady_Gradient_Spin_Echo_XY-16_Detection/op_sgse16_res.py @@ -0,0 +1,201 @@ +# -*- coding: iso-8859-1 -*- + +from numpy import * +from scipy.signal import * +from scipy.optimize import * +from os import path, rename + +def result(): + + measurement = MeasurementResult('Magnetization') + + suffix = '' # output file name's suffix and... + counter = 1 # counter for arrayed experiments + var_key = '' # variable parameter name + + # loop over the incoming results: + for timesignal in results: + if not isinstance(timesignal,ADC_Result): + continue + + # read experiment parameters: + pars = timesignal.get_description_dictionary() + + # get number of echoes: + num_echoes = pars['NECH'] + + # get phase list for data routung: + rec_phase = eval(pars['rec_phase']) + + # phase individual echoes in timesignal: + for i in range(num_echoes): + echo = timesignal.get_result_by_index(i) + echo.phase(rec_phase[i]) + + start = timesignal.index[i][0] + end = timesignal.index[i][1] + for k in range(timesignal.get_number_of_channels()): + if i%2 == 0: + timesignal.y[k][start:end+1] = echo.y[k] + else: + timesignal.y[k][start:end+1] = -echo.y[k] + + # provide timesignal to the display tab: + data['Current scan'] = timesignal + + # accumulate... + if not locals().get('accu'): + accu = Accumulation() + + # skip dummy scans, if any: + if pars['run'] < 0: continue + + # add up: + accu += timesignal + + # provide accumulation to the display tab: + data['Accumulation'] = accu + + # check how many scans are done: + if accu.n == pars['NS']: # accumulation is complete + + # get number of echoes: + num_echoes = pars['NECH'] + + # downsize accu to one point per echo: + echodecay = accu + 0 + echodecay.x = resize(echodecay.x, int(num_echoes)) + echodecay.y[0] = resize(echodecay.y[0], int(num_echoes)) + echodecay.y[1] = resize(echodecay.y[1], int(num_echoes)) + + # find where first echo exceeds noise: +# if not locals().get('noise'): +# last_echo = accu.get_accu_by_index(num_echoes-1) +# first_echo = accu.get_accu_by_index(0) +# noise = 3*std(last_echo.y[1][:]) +# samples = first_echo.y[0] > noise + + if not locals().get('noise'): + echo = accu.get_accu_by_index(0) + noise = 0.1*max(abs(echo.y[0])) + samples = abs(echo.y[0]) > noise + + # set echo times and intensities: + for i in range(num_echoes): + # get ith echo: + echo = accu.get_accu_by_index(i) + # set echo time: + echodecay.x[i] = i*2*pars['TAU'] + # set echo value: + echodecay.y[0][i] = sum(echo.y[0][samples]) # the sum of echo points that exceed noise... + echodecay.y[1][i] = sum(echo.y[1][samples]) + #echodecay.y[0][i] = sum(echo.y[0]) # or the sum of all echo points... + #echodecay.y[1][i] = sum(echo.y[1]) + #echodecay.y[0][i] = echo.y[0][echo.x.size/2] # or a middle echo point + #echodecay.y[1][i] = echo.y[1][echo.x.size/2] + + # compute a signal's phase: + phi0 = arctan2(echodecay.y[1][0], echodecay.y[0][0]) * 180 / pi + if not locals().get('ref'): ref = phi0 + print 'phi0 = ', phi0 + + # rotate signal to maximize Re (optional): + #echodecay.phase(-phi0) + + # provide echo decay to the display tab: + data['Echo Decay'] = echodecay + + # check whether it is an arrayed experiment: + var_key = pars.get('VAR_PAR') + if var_key: + # get variable parameter's value: + var_value = pars.get(var_key) + + # provide signal recorded with different var_value's to the display tab: + data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu + data['Echo Decay'+"/"+var_key+"=%e"%(var_value)] = echodecay + + # measure signal intensity vs. var_value: +# [amplitude, rate] = fitfunc(echodecay.x, echodecay.y[0]) +# print '%s%02g' % ('Amplitude = ', amplitude) +# print '%s%02g' % ('T2 [s] = ', 1./rate) + +# measurand[var_value] = amplitude + measurement[var_value] = sum(echodecay.y[0][:]) + + # provide measurement to the display tab: + data[measurement.get_title()] = measurement + + # save accu if required: + outfile = pars.get('OUTFILE') + if outfile: + datadir = pars.get('DATADIR') + + # write raw data in Simpson format: + filename = datadir+outfile+suffix+'.dat' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.dat') + accu.write_to_simpson(filename) + + # write raw data in Tecmag format: + # filename = datadir+outfile+'.tnt' + # accu.write_to_tecmag(filename, nrecords=20) + + # write parameters in a text file: + filename = datadir+outfile+suffix+'.par' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.par') + + fileobject = open(filename, 'w') + for key in sorted(pars.iterkeys()): + if key=='run': continue + if key=='rec_phase': continue + fileobject.write('%s%s%s'%(key,'=', pars[key])) + fileobject.write('\n') + fileobject.close() + + # reset accumulation: + del accu + + if var_key == 'D2': + # mono-exponential decay fit: + xdata = measurement.get_xdata() + ydata = measurement.get_ydata() + [amplitude, rate] = fitfunc(xdata, ydata) + print '%s%02g' % ('Amplitude = ', amplitude) + print '%s%02g' % ('T2 [s] = ', 1./rate) + + # update display for the fit: + measurement.y = func([amplitude, rate], xdata) + data[measurement.get_title()] = measurement + +# the fitting procedure: +def fitfunc(xdata, ydata): + + # initialize variable parameters: + try: + # solve Az = b: + A = array((ones(xdata.size/2), xdata[0:xdata.size/2])) + b = log(abs(ydata[0:xdata.size/2])) + z = linalg.lstsq(transpose(A), b) + amplitude = exp(z[0][0]) + rate = -z[0][1] + print amplitude, 1./rate + except: + amplitude = abs(ydata[0]) + rate = 1./(xdata[-1] - xdata[0]) + p0 = [amplitude, rate] + + # run least-squares optimization: + plsq = leastsq(residuals, p0, args=(xdata, ydata)) + + return plsq[0] # best-fit parameters + +def residuals(p, xdata, ydata): + return ydata - func(p, xdata) + +# here is the function to fit: +def func(p, xdata): + return p[0]*exp(-p[1]*xdata) + +pass \ No newline at end of file diff --git a/Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_exp.py b/Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_exp.py new file mode 100644 index 0000000..7f59f98 --- /dev/null +++ b/Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_exp.py @@ -0,0 +1,173 @@ +# -*- coding: iso-8859-1 -*- + +TXEnableDelay = 1e-6 +TXEnableValue = 0b0001 # TTL line blanking RF amplifier (bit 0) +TXPulseValue = 0b0010 # TTL line triggering RF pulses (bit 1) +ADCSensitivity = 1 # voltage span for ADC + +def experiment(): # Four-pulse STE sequence (Zeeman order) + + # set up acquisition parameters: + pars = {} + pars['P90'] = 4.2e-6 # 90-degree pulse length (s) + pars['SF'] = 338.7e6 # spectrometer frequency (Hz) + pars['O1'] = -60e3 # offset from SF (Hz) + pars['SW'] = 10e6 # spectral window (Hz) + pars['SI'] = 1*1024 # number of acquisition points + pars['NS'] = 16 # number of scans + pars['DS'] = 0 # number of dummy scans + pars['RD'] = .5 # delay between scans (s) + pars['D1'] = 30e-6 # delay after first pulse, or 'short tau' (s) + pars['D2'] = 100e-6 # delay after second pulse, or 'long tau' (s) + pars['D3'] = 20e-6 # refocusing pusle delay (s) + pars['PHA'] = 124 # receiver phase (degree) + pars['DATADIR'] = '/home/fprak/Students/' # data directory + pars['OUTFILE'] = None # output file name + + # specify a variable parameter (optional): + pars['VAR_PAR'] = 'D2' # variable parameter name (a string) + start = 30e-6 # starting value + stop = 2e-0 # end value + steps = 24 # number of values + log_scale = True # log scale flag + stag_range = False # staggered range flag + + # check parameters for safety: + if pars['PHA'] < 0: + pars['PHA'] = 360 + pars['PHA'] + + if pars['P90'] > 20e-6: + raise Exception("Pulse too long!!!") + + # check whether a variable parameter is named: + var_key = pars.get('VAR_PAR') + if var_key == 'P90' and (start > 20e-6 or stop > 20e-6): + raise Exception("Pulse too long!!!") + + if pars['NS']%8 != 0: + pars['NS'] = int(round(pars['NS'] / 8) + 1) * 8 + print 'Number of scans changed to ',pars['NS'],' due to phase cycling' + + # start the experiment: + if var_key: + # this is an arrayed experiment: + if log_scale: + array = log_range(start,stop,steps) + else: + array = lin_range(start,stop,steps) + + if stag_range: + array = staggered_range(array, size = 2) + + # estimate the experiment time: + if var_key == 'D1': + seconds = (sum(array)*2 + (pars['D2'] + pars['RD']) * steps) * (pars['NS'] + pars['DS']) + elif var_key == 'D2': + seconds = (sum(array) + (pars['D1']*2 + pars['RD']) * steps) * (pars['NS'] + pars['DS']) + elif var_key == 'RD': + seconds = (sum(array) + (pars['D1']*2 + pars['D2']) * steps) * (pars['NS'] + pars['DS']) + else: + seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * steps * (pars['NS']+ pars['DS']) + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for a variable parameter: + for index, pars[var_key] in enumerate(array): + print 'Arrayed experiment for '+var_key+': run = '+str(index+1)+\ + ' out of '+str(array.size)+': value = '+str(pars[var_key]) + # loop for accumulation: + for run in xrange(pars['NS']+pars['DS']): + yield zeeman4pulses_experiment(pars, run) + synchronize() + else: + # estimate the experiment time: + seconds = (pars['D1']*2 + pars['D2'] + pars['RD']) * (pars['NS']+ pars['DS']) + m, s = divmod(seconds, 60) + h, m = divmod(m, 60) + print '%s%02d:%02d:%02d' % ('Experiment time estimated: ', h, m, s) + + # loop for accumulation: + for run in xrange(pars['NS']+pars['DS']): + yield zeeman4pulses_experiment(pars, run) + +# the pulse program: + +def zeeman4pulses_experiment(pars, run): + e=Experiment() + + dummy_scans = pars.get('DS') + if dummy_scans: + run -= dummy_scans + + pars['PROG'] = 'zeeman4pulses_experiment' + + # ok 8-step phase cycle (1-21 modifided to deal with T1-recovery and extended for Re/Im imbalance) + pars['PH1'] = [0, 270, 0, 270, 180, 90, 180, 90] # 1st (90-degree) pulse + pars['PH3'] = [0, 90, 0, 90, 0, 90, 0, 90] # 2nd (90-degree) pulse + pars['PH4'] = [0, 0, 180, 180, 270, 270, 90, 90] # 3rd (90-degree) pulse + pars['PH5'] = [90, 90, 90, 90, 180, 180, 180, 180] # 3rd (90-degree) pulse + pars['PH2'] = [0, 180, 180, 0, 90, 270, 270, 90] # receiver + + + # read in variables: + P90 = pars['P90'] + SF = pars['SF'] + O1 = pars['O1'] + RD = pars['RD'] + D1 = pars['D1'] + D2 = pars['D2'] + D3 = pars['D3'] + PH1 = pars['PH1'][run%len(pars['PH1'])] + PH3 = pars['PH3'][run%len(pars['PH3'])] + PH4 = pars['PH4'][run%len(pars['PH4'])] + PH5 = pars['PH5'][run%len(pars['PH5'])] + PH2 = pars['PH2'][run%len(pars['PH2'])] + PHA = pars['PHA'] + + if (run/8)%2 != 0: + PH5 += 180 + + # set sampling parameters: + SI = pars['SI'] + SW = pars['SW'] + while SW <= 10e6 and SI < 256*1024: + SI *= 2 + SW *= 2 + + # run the pulse program: + e.wait(RD) # relaxation delay between scans + e.set_frequency(SF+O1, phase=PH1) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(D1-P90-TXEnableDelay) # short tau + + e.set_phase(PH3) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(D2-P90-TXEnableDelay) # long tau + + e.set_phase(PH4) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(D3-P90-TXEnableDelay) # echo delay + + e.set_phase(PH5) + e.ttl_pulse(TXEnableDelay, value=TXEnableValue) + e.ttl_pulse(P90, value=TXEnableValue|TXPulseValue) # 90-degree pulse + + e.wait(TXEnableDelay) + e.set_phase(PHA) + e.wait(D3+D1-P90/2-TXEnableDelay) # echo delay + e.record(SI, SW, sensitivity=ADCSensitivity) # acquisition + + # write experiment parameters: + for key in pars.keys(): + e.set_description(key, pars[key]) # acquisition parameters + e.set_description('run', run) # current scan + e.set_description('rec_phase', -PH2) # current receiver phase + + return e \ No newline at end of file diff --git a/Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_res.py b/Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_res.py new file mode 100644 index 0000000..3f5cb26 --- /dev/null +++ b/Scripts/Zeeman_Order_Four_Pulses/op_zeeman4pulses_res.py @@ -0,0 +1,210 @@ +# -*- coding: iso-8859-1 -*- + +from numpy import * +from scipy.signal import * +from scipy.optimize import * +from os import path, rename + +def result(): + + measurement = MeasurementResult('Magnetization') + + measurement_range = [0.0, 10e-6] + measurement_ranging = False + + suffix = '' # output file name's suffix and... + counter = 1 # counter for arrayed experiments + var_key = '' # variable parameter name + + # loop over the incoming results: + for timesignal in results: + if not isinstance(timesignal,ADC_Result): + continue + + # read experiment parameters: + pars = timesignal.get_description_dictionary() + + # ---------------- digital filter ------------------ + + # get actual sampling rate of timesignal: + sampling_rate = timesignal.get_sampling_rate() + + # get user-defined spectrum width: + spec_width = pars['SW'] + + # specify cutoff frequency, in relative units: + cutoff = spec_width / sampling_rate + + if cutoff < 1: # otherwise no filter applied + + # number of filter's coefficients: + numtaps = 29 + + # use firwin to create a lowpass FIR filter: + fir_coeff = firwin(numtaps, cutoff) + + # downsize x according to user-defined spectral window: + skip = int(sampling_rate / spec_width) + timesignal.x = timesignal.x[::skip] + + for i in range(2): + # apply the filter to ith channel: + timesignal.y[i] = lfilter(fir_coeff, 1.0, timesignal.y[i]) + + # zeroize first N-1 "corrupted" samples: + timesignal.y[i][:numtaps-1] = 0.0 + + # circular left shift of y: + timesignal.y[i] = roll(timesignal.y[i], -(numtaps-1)) + + # downsize y to user-defined number of samples (SI): + timesignal.y[i] = timesignal.y[i][::skip] + + # update the sampling_rate attribute of the signal's: + timesignal.set_sampling_rate(spec_width) + + # ---------------------------------------------------- + + # rotate timesignal according to current receiver's phase: + timesignal.phase(pars['rec_phase']) + + # provide timesignal to the display tab: + data['Current scan'] = timesignal + + # accumulate... + if not locals().get('accu'): + accu = Accumulation() + + # skip dummy scans, if any: + if pars['run'] < 0: continue + + # add up: + accu += timesignal + + # provide accumulation to the display tab: + data['Accumulation'] = accu + + # check how many scans are done: + if accu.n == pars['NS']: # accumulation is complete + + # make a copy: + echo = accu + 0 + + # compute the initial phase of FID: + phi0 = arctan2(echo.y[1][0], echo.y[0][0]) * 180 / pi + if not 'ref' in locals(): ref = phi0 + print 'phi0 = ', phi0 + + # rotate FID to maximize y[0][0]: + #echo.phase(-phi0) + + # do FFT: + echo.exp_window(line_broadening=10) + spectrum = echo.fft(samples=2*pars['SI']) + + # try zero-order phase correction: + spectrum.phase(-phi0) + + # provide spectrum to the display tab: + data['Spectrum'] = spectrum + + # check whether it is an arrayed experiment: + var_key = pars.get('VAR_PAR') + if var_key: + # get variable parameter's value: + var_value = pars.get(var_key) + + # provide signal recorded with this var_value to the display tab: + data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu + + # measure signal intensity vs. var_value: + if measurement_ranging == True: + [start, stop] = accu.get_sampling_rate() * array(measurement_range) + measurement[var_value] = sum(accu.y[0][int(start):int(stop)]) + + else: + measurement[var_value] = sum(accu.y[0][0:31]) + + # provide measurement to the display tab: + data[measurement.get_title()] = measurement + + # update the file name suffix: + suffix = '_' + str(counter) + counter += 1 + + # save accu if required: + outfile = pars.get('OUTFILE') + if outfile: + datadir = pars.get('DATADIR') + + # write raw data in Simpson format: + filename = datadir+outfile+suffix+'.dat' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.dat') + accu.write_to_simpson(filename) + + # write raw data in Tecmag format: + # filename = datadir+outfile+'.tnt' + # accu.write_to_tecmag(filename, nrecords=20) + + # write parameters in a text file: + filename = datadir+outfile+suffix+'.par' + if path.exists(filename): + rename(filename, datadir+'~'+outfile+suffix+'.par') + + fileobject = open(filename, 'w') + for key in sorted(pars.iterkeys()): + if key=='run': continue + if key=='rec_phase': continue + fileobject.write('%s%s%s'%(key,'=', pars[key])) + fileobject.write('\n') + fileobject.close() + + # reset accumulation: + del accu + + + if var_key == 'D2': + # KWW fit: + xdata = measurement.get_xdata() + ydata = measurement.get_ydata() + [amplitude, rate, beta] = fitfunc(xdata, ydata) + print '%s%02g' % ('Amplitude = ', amplitude) + print '%s%02g' % ('T2 [s] = ', 1./rate) + print '%s%02g' % ('Beta = ', beta) + + # update display for the fit: + measurement.y = func([amplitude, rate, beta], xdata) + data[measurement.get_title()] = measurement + +# the fitting procedure: +def fitfunc(xdata, ydata): + + # initialize variable parameters: + try: + # solve Az = b: + A = array((ones(xdata.size/2), xdata[0:xdata.size/2])) + b = log(abs(ydata[0:xdata.size/2])) + z = linalg.lstsq(transpose(A), b) + amplitude = exp(z[0][0]) + rate = -z[0][1] + except: + amplitude = abs(ydata[0]) + rate = 1./(xdata[-1] - xdata[0]) + beta = 1 + p0 = [amplitude, rate, beta] + + # run least-squares optimization: + plsq = leastsq(residuals, p0, args=(xdata, ydata)) + + return plsq[0] # best-fit parameters + +def residuals(p, xdata, ydata): + return ydata - func(p, xdata) + +# here is the function to fit: +def func(p, xdata): + return p[0]*exp(-(p[1]*xdata)**p[2]) + + +pass \ No newline at end of file