pade works

phase-encoding-method-of-moments-Pade.py → pade-method-of-moments.py

@@ -17,8 +17,8 @@ U_points = voxel_num * 3000
 
 noise_ratio = 0.0*(1.2e-16) #(1.2e-16) #1e8
 
-total_periods=10#DO NOT CHANGE to fit T2s_scale autoscale
-rf_samples_per_period = 4
+total_periods=5#10#DO NOT CHANGE to fit T2s_scale autoscale
+rf_samples_per_period = 2
 
 # max polynomial order equals  rf_samples_per_period * total_periods 
 
@@ -30,7 +30,8 @@ T2s_scale = total_periods*period/rf_samples_per_period/voxel_num*8 #ms # need to
 T2s_min = T2s_scale/20.0 
 #ms
 total_steps = rf_samples_per_period*total_periods
-if total_steps % 2 != 0: total_steps -= 1
+n_mu = total_steps
+#if total_steps % 2 != 0: total_steps -= 1
 time_steps = np.array(mp.linspace((0), (period*total_periods), total_steps))
 tmp = [np.random.rand()*0.0+(1.0) for n in range(voxel_num)]
 voxel_amplitudes = np.array(tmp)
@@ -45,9 +46,9 @@ if voxel_num == 5:
     voxel_amplitudes = voxel_all[:voxel_num]
     voxel_T2s_decay = voxel_all[voxel_num:]*T2s_scale
 
-# a_i = np.array([(0.3),(0.1),(0.15),(0.1)])
-# d_i = np.array([(0.7),(0.9),(0.2),(0.67)])
-# voxel_num = len(a_i)
+a_i = np.array([(0.3),(0.1),(0.15),(0.1),(0.2)])
+d_i = np.array([(0.7),(0.9),(0.2),(0.67),(0.41)])
+voxel_num = len(a_i)
 
 
 voxel_phases = np.array(np.linspace(0,phase_range, voxel_num))
@@ -63,23 +64,34 @@ 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)
+    v_sum = 0
     for t in range(len(time)):
         for i in range(voxel_num):
-            # mag_sin[t] += a_i[i]*(
-            #     (d_i[i] + np.random.rand()*noise_ratio)**time[t]
-            # )
-            amp = voxel_amplitudes[i] * (
-                mp.exp(-time[t]/voxel_T2s_decay[i])
-                ) + ( 0.0 
-                     #     + np.random.rand()*noise_ratio
-                        )
+            mag_sin[t] += a_i[i]*(
+                (d_i[i] + np.random.rand()*noise_ratio)**t
+            )
+            # 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]* np.sin(voxel_phases[i] + phase_init)))
+            #     #print("d_{:d}".format(i),float( np.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
+            # mag_sin[t] += amp * np.sin((voxel_phases[i] + phase_init))
+            # mag_cos[t] += amp * np.cos((voxel_phases[i] + phase_init))
             # if t == 0:
-                #print("a_{:d}".format(i),float(voxel_amplitudes[i]* np.sin(voxel_phases[i] + phase_init)))
-                #print("d_{:d}".format(i),float( np.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
-            mag_sin[t] += amp * np.sin((voxel_phases[i] + phase_init))
-            mag_cos[t] += amp * np.cos((voxel_phases[i] + phase_init))
+            #     v = voxel_amplitudes[i] * (
+            #     mp.exp(-(time[1])/voxel_T2s_decay[i])
+            #     )* np.sin((voxel_phases[i] + phase_init))
+            #     v_sum += v
+            #     print(float(v))
+            #     #sign+="{:4.3g}, ".format(float( mp.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
+
         mag_sin[t] *= (1 + (np.random.rand()-0.5)*noise_ratio)
+        # if t == 0:
+        #     print(v_sum, mag_sin[t])
         #mag_cos[t] += np.random.rand()*noise_ratio
     return mag_sin, mag_cos
 
@@ -88,9 +100,51 @@ mag_sin_0, mag_cos_0 = gen_rf_signal(time_steps)
 mag_sin = np.copy(mag_sin_0)
 mag_cos = np.copy(mag_cos_0)
 
+# mu = [np.random.rand() for i in range(n_mu)]
+mu = mag_sin
+#print(mu)
+
+P = np.zeros((n_mu+1, n_mu+1),dtype=mpf)
+P[0,0] = 1
+for i in range(n_mu):
+    P[i,1] = (-1)**i * mu[i]
+
+for j in range(2,n_mu+1):
+    for i in range(n_mu-1,-1,-1):
+        if i+j > n_mu: continue
+        # print("i =",i, "  j =",j, "   P = ",P[i,j])
+        P[i,j] = P[0,j-1]*P[i+1,j-2] - P[0,j-2]*P[i+1,j-1]
+
+# for i in range (n_mu+1):
+#     line = [float("%5.3g"%(x)) for x in P[i,:]]
+#     print(line)
+#print(P)
 
-print("N = ",len(mag_sin))
-print("Noise = ", noise_ratio )
+alpha = np.zeros(n_mu,dtype=mpf)
+for k in range(1,n_mu):
+    alpha[k] = (P[0,k+1])/((P[0,k])*(P[0,k-1]))
+# print("alpha:",len(alpha))
+# line = [float("%5.5g"%(x)) for x in alpha]
+# print(line)
+
+
+M=np.zeros((voxel_num,voxel_num), dtype=mpf)
+for i in range(voxel_num):
+    for j in range(voxel_num):
+        if i == j:
+            M[i,j] = alpha[i*2]+alpha[i*2 +1]
+        if abs(i-j) == 1:
+            m = max(i,j)
+            M[i,j] = -np.sqrt( alpha[2*m-1] * alpha[2*m] )
+M1 = mp.matrix(M.tolist())
+#print(M1)
+
+spectra = mp.eig(M1)[0]
+line = [float("%5.5g"%(x)) for x in spectra]
+print("eig(M1) =",line)
+
+# print("N = ",len(mag_sin))
+# print("Noise = ", noise_ratio )
 # print("c_n =",["%0.8f"%(x) for x in mag_sin])
 
 sign = ""
@@ -108,12 +162,12 @@ import time
 time1 = time.time()
 
 time2 = time.time()
-print( 'It took %0.5f s' % ((time2-time1)))
+#print( 'It took %0.5f s' % ((time2-time1)))
 
 sign ="" 
 for i in range(voxel_num):
     sign+="{:4.3g}, ".format(float( mp.exp(-(period/rf_samples_per_period)/voxel_T2s_decay[i]) ))
-#print(sign)
+print(d_i)
 # x = np.linspace(0,1, U_points)
 # mag_x = np.linspace(0,1, len(mag_sin))
 # import matplotlib as matplotlib