method of moment works

phase-encoding-method-of-moments.py

@@ -15,34 +15,34 @@ U_points = voxel_num * 1000
 # noise_ratio = mpf(0.0) #1e8
 
 total_periods = 10
-rf_samples_per_period = 5
+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(15.6/1000/1000) #ms
+period = mpf(1/(total_periods*rf_samples_per_period)) #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
 #print(period)
 #ms
-time_steps = np.array(mp.linspace(mpf(0), mpf(period*total_periods), rf_samples_per_period*total_periods))
+time_steps = np.array(mp.linspace(mpf(0), mpf(rf_samples_per_period*total_periods), rf_samples_per_period*total_periods+1))
 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_all = np.append(voxel_amplitudes,voxel_T2s_decay/T2s_scale)
-if voxel_num == 5:
-    voxel_all = np.array([mpf(0.822628),mpf(0.691376),mpf(0.282906),mpf(0.226013),mpf(0.90703),mpf(0.144985),mpf(0.328563),mpf(0.440353),mpf(0.662462),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.6)])
+voxel_num = len(a_i)
+
 
 voxel_phases = np.array(mp.linspace(0,phase_range, voxel_num))
-if len(voxel_amplitudes) != len(voxel_phases):
-    print("ERROR! Size of amplitude and phase arrays do not match!")
-    raise
+# if len(voxel_amplitudes) != len(voxel_phases):
+#     print("ERROR! Size of amplitude and phase arrays do not match!")
+#     raise
 
+ampl = []
 def gen_rf_signal(time):
     '''Generates demodulated signal at radio frequence using voxels
     amplitudes, T2s decays, and phases.
@@ -50,20 +50,29 @@ def gen_rf_signal(time):
     '''
     tmp = [mpf(0.0) for n in range(len(time))]
     mag_sin = np.array(tmp)
-    mag_cos = np.array(tmp)
+    mag_cos = np.array(tmp)    
     for t in range(len(time)):
+        #print("time",float(time[t]))
         for i in range(voxel_num):
-            amp = voxel_amplitudes[i] * (
-                mp.exp(-time[t]/voxel_T2s_decay[i])
-                ) + ( 0.0 
-                          # + np.random.rand()*noise_ratio
-                        )
-            mag_sin[t] += amp * mp.sin(
-                voxel_phases[i] + phase_init
-                )
-            mag_cos[t] += amp * mp.cos(
-                voxel_phases[i] + phase_init
-                )
+            mag_sin[t] += a_i[i]*(d_i[i]**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
+            #     )
     return mag_sin, mag_cos
 
 def factorial(n):
@@ -76,7 +85,7 @@ def shiftedLegendre(n):
     coeffs = []
     for k in range(n+1):
         val = mpf(-1)**n * binom(mpf(n),mpf(k)) * binom(n+k,k) * (-1)**k
-        coeffs.insert(0,val)
+        coeffs.insert(0,val*mp.sqrt(2*n+1))
     return np.poly1d(coeffs)
 
 def K ( i, j):
@@ -98,13 +107,15 @@ def GetLambda(mag_rf):
     M_cutoff = len(mag_rf)
     all_lambda = []
     for i in range(M_cutoff):
+#        print("M = ", i)
         lambd = mpf(0)
         for j in range(i+1):
-            lambd += K(i,j)*mag_rf[i]
+            lambd += K(i,j)*mag_rf[j]
+#            print("K({:d},{:d}) =".format(i,j), float(K(i,j)))
         all_lambda.append(lambd)
     # tmp = [mpf(0.0) for n in range(M_cutoff)]
     # all_lambda = np.array(tmp)
-    # all_lambda[29] = mpf(1.0)
+    # all_lambda[10] = mpf(1.0)
     return all_lambda
 
 
@@ -114,12 +125,12 @@ sign = ""
 for i in range(voxel_num):
     if i%5 == 0:
         sign+="\n"
-    sign = sign + '{:3.2g}'.format(float(voxel_amplitudes[i] * mp.sin(
-        voxel_phases[i] + phase_init
-        )))+"/"+'{:3.2g}'.format(float(voxel_T2s_decay[i]))+", "
-
+    sign = sign + '{:3.2g}'.format(float(a_i[i]))+"/"+'{:3.2g}'.format(float(d_i[i]))+", "
 
-plt.plot(mag_sin, ls='-', marker='o')
+#    print(mp.exp(-1.0/voxel_T2s_decay[i]))
+               
+ 
+plt.plot(mag_sin, ls=' ', marker='o')
 plt.title("Signal to restore amp/decay_T:"+sign)
 plt.savefig("signal.pdf")
 plt.clf()
@@ -127,7 +138,9 @@ plt.clf()
 
 L = [] # Shifted Legendre polinomials
 for i in range(len(mag_sin)):
-    polyL = shiftedLegendre(i)        
+    polyL = shiftedLegendre(i)
+    # print("i=",i,"  L_i:")
+    # print(polyL)
     L += [polyL]
 
 x = np.linspace(0,1, U_points)
@@ -147,6 +160,7 @@ U = GetU(lambdas)
 x = np.linspace(0,1, U_points)
 mag_x = np.linspace(0,1, len(mag_sin))
 plt.plot(x,U)
+plt.title(r"$\sum^M \lambda_i L_i(x)$", y=1.02)
 plt.savefig("plt.pdf")
 plt.clf()