|
|
@@ -5,36 +5,44 @@ import numpy as np
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
from mpmath import mp, mpf
|
|
|
-mp.dps = 1000
|
|
|
+mp.dps = 200
|
|
|
|
|
|
-voxel_num = 2
|
|
|
+voxel_num = 5
|
|
|
phase_range = mp.pi/2
|
|
|
-phase_init = mp.pi/20 #mpf(0.0)
|
|
|
+phase_init = mp.pi/4 #mpf(0.0)
|
|
|
U_points = voxel_num * 1000
|
|
|
|
|
|
-# noise_ratio = mpf(0.0) #1e8
|
|
|
+noise_ratio = mpf(0.0*1e-38) #1e8
|
|
|
|
|
|
total_periods = 10
|
|
|
rf_samples_per_period = 10
|
|
|
# max polynomial order equals rf_samples_per_period * total_periods
|
|
|
|
|
|
# B0=1.5T freq=64Mhz, period = 15.6 ns
|
|
|
-period = mpf(1/(total_periods*rf_samples_per_period)) #ms
|
|
|
+period = mpf(10) #ms
|
|
|
omega = 2.0*mp.pi/period
|
|
|
#T2s_scale = 0.01 #ms # need to be 10ms
|
|
|
-T2s_scale = total_periods*period #ms # need to be 10ms
|
|
|
-T2s_min = T2s_scale/1000.0
|
|
|
+T2s_scale = total_periods*period/15 #ms # need to be 10ms
|
|
|
+T2s_min = T2s_scale/10.0
|
|
|
#print(period)
|
|
|
#ms
|
|
|
-time_steps = np.array(mp.linspace(mpf(0), mpf(rf_samples_per_period*total_periods), rf_samples_per_period*total_periods+1))
|
|
|
+time_steps = np.array(mp.linspace(mpf(0), mpf(period*total_periods), rf_samples_per_period*total_periods))
|
|
|
tmp = [mp.rand() for n in range(voxel_num)]
|
|
|
voxel_amplitudes = np.array(tmp)
|
|
|
tmp = [mp.rand() for n in range(voxel_num)]
|
|
|
-voxel_T2s_decay = np.array(tmp)*(T2s_scale-T2s_min) + T2s_min
|
|
|
+voxel_T2s_decay = np.array(tmp)*(T2s_scale-2*T2s_min) + T2s_min
|
|
|
+print(voxel_T2s_decay)
|
|
|
voxel_all = np.append(voxel_amplitudes,voxel_T2s_decay/T2s_scale)
|
|
|
-a_i = np.array([mpf(0.3),mpf(0.1),mpf(0.15),mpf(0.1)])
|
|
|
-d_i = np.array([mpf(0.7),mpf(0.9),mpf(0.2),mpf(0.6)])
|
|
|
-voxel_num = len(a_i)
|
|
|
+if voxel_num == 5:
|
|
|
+# voxel_all = np.array([mpf(0.2),mpf(0.6),mpf(0.02),mpf(0.1)])
|
|
|
+ voxel_all = np.array([mpf(0.822628),mpf(0.691376),mpf(0.282906),mpf(0.226013),mpf(0.90703),mpf(0.144985),mpf(0.228563),mpf(0.340353),mpf(0.462462),mpf(0.720518)])
|
|
|
+ #voxel_all = np.array([mpf(0.592606),mpf(0.135168),mpf(0.365712),mpf(0.667536),mpf(0.437378),mpf(0.918822),mpf(0.943879),mpf(0.590338),mpf(0.685997),mpf(0.658839)])
|
|
|
+ voxel_amplitudes = voxel_all[:voxel_num]
|
|
|
+ voxel_T2s_decay = voxel_all[voxel_num:]*T2s_scale
|
|
|
+
|
|
|
+# a_i = np.array([mpf(0.3),mpf(0.1),mpf(0.15),mpf(0.1)])
|
|
|
+# d_i = np.array([mpf(0.7),mpf(0.9),mpf(0.2),mpf(0.67)])
|
|
|
+# voxel_num = len(a_i)
|
|
|
|
|
|
|
|
|
voxel_phases = np.array(mp.linspace(0,phase_range, voxel_num))
|
|
|
@@ -54,25 +62,26 @@ def gen_rf_signal(time):
|
|
|
for t in range(len(time)):
|
|
|
#print("time",float(time[t]))
|
|
|
for i in range(voxel_num):
|
|
|
- mag_sin[t] += a_i[i]*(d_i[i]**time[t])
|
|
|
+ # mag_sin[t] += a_i[i]*(
|
|
|
+ # (d_i[i] + np.random.rand()*noise_ratio)**time[t]
|
|
|
+ # )
|
|
|
# print("a_{:d} =".format(i),float(a_i[i]),", ",
|
|
|
- # "d_{:d} =".format(i),float(d_i[i]))
|
|
|
-
|
|
|
- # amp = voxel_amplitudes[i] * (
|
|
|
- # mp.exp(-time[t]/voxel_T2s_decay[i])
|
|
|
- # ) + ( 0.0
|
|
|
- # # + np.random.rand()*noise_ratio
|
|
|
- # )
|
|
|
- # print("a_{:d}".format(i),float(voxel_amplitudes[i]* mp.sin(
|
|
|
- # voxel_phases[i] + phase_init
|
|
|
- # )))
|
|
|
- # print("d_{:d}".format(i),float( mp.exp(-1.0/voxel_T2s_decay[i]) ))
|
|
|
- # mag_sin[t] += amp * mp.sin(
|
|
|
- # voxel_phases[i] + phase_init
|
|
|
- # )
|
|
|
- # mag_cos[t] += amp * mp.cos(
|
|
|
- # voxel_phases[i] + phase_init
|
|
|
- # )
|
|
|
+ # "d_{:d} =".format(i),float(d_i[i]),"+", np.random.rand()*noise_ratio)
|
|
|
+
|
|
|
+ amp = voxel_amplitudes[i] * (
|
|
|
+ mp.exp(-time[t]/voxel_T2s_decay[i])
|
|
|
+ ) + ( 0.0
|
|
|
+ # + np.random.rand()*noise_ratio
|
|
|
+ )
|
|
|
+ if t == 0:
|
|
|
+ #print("a_{:d}".format(i),float(voxel_amplitudes[i]* mp.sin(voxel_phases[i] + phase_init)))
|
|
|
+ print("d_{:d}".format(i),float( mp.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
|
|
|
+ mag_sin[t] += amp * mp.sin(
|
|
|
+ voxel_phases[i] + phase_init
|
|
|
+ )
|
|
|
+ mag_cos[t] += amp * mp.cos(
|
|
|
+ voxel_phases[i] + phase_init
|
|
|
+ )
|
|
|
return mag_sin, mag_cos
|
|
|
|
|
|
def factorial(n):
|
|
|
@@ -99,6 +108,7 @@ def GetU (lambdas):
|
|
|
tmp = [mpf(0.0) for n in range(U_points)]
|
|
|
U = np.array(tmp)
|
|
|
for i in range (len(lambdas)):
|
|
|
+ if i%10 == 0: print(i)
|
|
|
polyL = L[i] #shiftedLegendre(i)
|
|
|
U += lambdas[i]*polyL(x)
|
|
|
return U
|
|
|
@@ -122,12 +132,12 @@ def GetLambda(mag_rf):
|
|
|
mag_sin, mag_cos = gen_rf_signal(time_steps)
|
|
|
|
|
|
sign = ""
|
|
|
-for i in range(voxel_num):
|
|
|
- if i%5 == 0:
|
|
|
- sign+="\n"
|
|
|
- sign = sign + '{:3.2g}'.format(float(a_i[i]))+"/"+'{:3.2g}'.format(float(d_i[i]))+", "
|
|
|
+# for i in range(voxel_num):
|
|
|
+# if i%5 == 0:
|
|
|
+# sign+="\n"
|
|
|
+# sign = sign + '{:3.2g}'.format(float(a_i[i]))+"/"+'{:3.2g}'.format(float(d_i[i]))+", "
|
|
|
|
|
|
-# print(mp.exp(-1.0/voxel_T2s_decay[i]))
|
|
|
+# # print(mp.exp(-1.0/voxel_T2s_decay[i]))
|
|
|
|
|
|
|
|
|
plt.plot(mag_sin, ls=' ', marker='o')
|
