Index: trunk/src/src/RICS_paper_compare_noise.py =================================================================== diff --git a/trunk/src/src/RICS_paper_compare_noise.py b/trunk/src/src/RICS_paper_compare_noise.py --- a/trunk/src/src/RICS_paper_compare_noise.py (revision 25) +++ b/trunk/src/src/RICS_paper_compare_noise.py (revision 26) @@ -11,6 +11,7 @@ import time import matplotlib.pyplot as plt from scipy.optimize import root +from scipy.optimize import newton_krylov from scipy.stats import nanmean from y_hysell96 import * @@ -53,38 +54,28 @@ ## Rotation matrix from antenna coord to magnetic coord (East North) -> 45 degrees theta_rad = np.radians(theta) # trig functions take radians as argument val1 = float( np.cos(theta_rad) ) -val2 = float( np.sin(theta_rad) ) -val3 = float( -1*np.sin(theta_rad)) +val2 = float( -1*np.sin(theta_rad)) +val3 = float( np.sin(theta_rad) ) val4 = float( np.cos(theta_rad) ) # Rotation matrix from antenna coord to magnetic coord (East North) -R = np.array( [[val1, val3], [val2, val4]] ) +R = np.array( [[val1, val2], [val3, val4]] ) # Rotate antenna positions to magnetic coord. AR = np.dot(R.T, antpos) -plt.plot(AR[0], 0*AR[1], 'yo') -plt.title("Jicamarca Antenna Positions") +plt.plot(AR[0], 0*AR[1], 'yo') # antenna positions only considering magnetic East-West component (x, 0) +plt.title("Jicamarca Antenna Positions") plt.xlabel("rx") plt.ylabel("ry") plt.grid(True) plt.plot(rx, 0*rx, 'g-') plt.draw() -print antpos -print "\n" -print AR # Only take the East component r = AR[0,:] r.sort() - -# ---NEW--- -r1 = AR[1,:] # longitudes / distances -d = 40 # separation of antennas in meters -theta_w = 15 # angle between antenna and point of observation in degrees -# arg = k*(r + d*np.cos(theta_w)) # this is the exponent argument in the corr integral -# -------- # Truth model (high and low resolution) @@ -123,7 +114,7 @@ fact = np.zeros(shape=(Nt,1)) factr = np.zeros(shape=(Nr,1)) -# Constructing Gaussian model +# Constructing Gaussian model - brightness function (?) for i in range(0, a.size): temp = (-(thetat-mu[i])**2/(sig[i]**2)) tempr = (-(thetar-mu[i])**2/(sig[i]**2)) @@ -170,11 +161,11 @@ # Plot Gaussian pulse(s) plt.figure(2) -plt.plot(thetat, fact, 'r--') +plt.plot(thetat, fact, 'bo') plt.plot(thetar, factr, 'ro') -plt.title("Truth models") +plt.title("Image and Reconstruction Domain") plt.ylabel("Amplitude") -plt.legend(('High res','Low res'), loc='upper right') +plt.legend(('High res (img domain)','Low res (reconstruction domain)'), loc='upper right') plt.draw() @@ -199,16 +190,67 @@ SNR = 10**(SNRdB/10.0); for Ni in range(0, NN): + + # -------------------------------------- Voltages estimation---------------------------------------- + v = np.zeros(shape=(r.size,1), dtype=object); # voltages matrix + R = np.zeros(shape=(r.size, r.size), dtype='complex128'); # corr matrix + mu1, sigma1 = 0, 1 # mean and standard deviation + num = 1 + + for t in range (0, 10): + num = t+1 + print num + n1 = np.random.normal(mu1, sigma1, [Nt,1]) # random variable + n2 = np.random.normal(mu1, sigma1, [Nt,1]) # random variable + + rho = (n1 + n2*1j)*np.sqrt(fact/2) # rho.size is Nt = 16834 + # fact is size (Nt, 1) -> column vector + thetat = thetat.reshape(Nt,1) # thetat is now (Nt, 1) + + for i in range(0, r.size): + v[i] = sum((rho)*np.exp(-1j*k*((r[i])*np.cos(thetat)))); # v is sum ( (Nt, 1)) = 1 value + + # Computing correlations + for i1 in range(0, r.size): + for i2 in range(0,r.size): + R[i1,i2] = R[i1,i2] + (v[i1][0]*v[i2][0].conjugate())/ (np.sqrt(np.abs(v[i1][0])**2 * np.abs(v[i2][0])**2)); + + R = R / num + print "R from voltages" + print R + + plt.figure(3); + plt.pcolor(abs(R)); + plt.colorbar(); + plt.title("Antenna Correlations") + plt.draw() + plt.show() + + #-------------------------------------------------------------------------------------------------- + # Calculate cross-correlation matrix (Fourier components of image) # This is an inefficient way to do this. - + + # MATLAB CORR::: R(i1,i2) = fact.' * exp(j*k*(r(i1)-r(i2))*sin(thetat)); + R = np.zeros(shape=(r.size, r.size), dtype=object); - + + thetat = thetat.reshape(Nt,) for i1 in range(0, r.size): for i2 in range(0,r.size): - R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat)))) - R[i1,i2] = sum(R[i1,i2]) - + R[i1,i2] = sum(np.dot(fact.T, np.exp(1j*k*((r[i1]-r[i2])*np.sin(thetat))))) + + print "R w/o voltages\n", R +# plt.figure(3); +# plt.pcolor(abs(R)); +# plt.colorbar(); +# plt.title("Antenna Correlations") +# plt.xlabel("Antenna") +# plt.ylabel("Antenna") +# plt.draw() +# plt.show() + + # Add uncertainty # This is an ad-hoc way of adding "noise". It models some combination of # receiver noise and finite integration times. We could use a more @@ -225,9 +267,9 @@ Unz = U + nz; Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R - plt.figure(3); - plt.pcolor(abs(R)); - plt.colorbar(); +# plt.figure(3); +# plt.pcolor(abs(R)); +# plt.colorbar(); #------------------------------------------------------------------------------------------------- # Fourier Inversion @@ -236,9 +278,8 @@ for i in range(0, thetar.size): th = thetar[i]; - w = np.exp(1j*k*np.dot(r,np.sin(th))); - temp = np.dot(w.T.conj(),U) - f_fourier[i] = np.dot(temp, w); + w = np.exp(1j*k*np.dot(r,np.sin(th))); + f_fourier[i] = np.dot(np.dot(w.T.conj(),U) , w); f_fourier = f_fourier.real; # get rid of numerical imaginary noise @@ -281,7 +322,7 @@ # Build H for i1 in range(0, r.size): for i2 in range(i1+1, r.size): - idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward 1->27 + idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward (1->27) idx1 = 2*idx; idx2 = 2*idx+1; Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj(); @@ -308,7 +349,8 @@ tic_maxEnt = time.time(); # start time maxEnt - lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14); + lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14); # krylov +# lambda1 = newton_krylov(myfun,lambda0, method='lgmres', x_tol=1e-14) toc_maxEnt = time.time() elapsed_time_maxent = toc_maxEnt - tic_maxEnt; @@ -319,11 +361,11 @@ f_maxent = modelf(lambda1, Hr, F); #------ what is this needed for?------- - ystar = myfun(lambda1); - Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G)); - ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda); - es = np.dot(Hr, f_maxent) - gnz; - chi2 = np.sum((es/sigma)**2); +# ystar = myfun(lambda1); +# Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G)); +# ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda); +# es = np.dot(Hr, f_maxent) - gnz; +# chi2 = np.sum((es/sigma)**2); #-------------------------------------- @@ -354,12 +396,13 @@ print "Psi ", Psi.shape print "s ", s.shape print "f_CS ", f_cs.shape + - # Plot - plt.figure(4) - plt.plot(thetar,f_cs,'r.-'); - plt.plot(thetat,fact,'k-'); - plt.title('Reconstruction with Compressed Sensing') +# # Plot +# plt.figure(4) +# plt.plot(thetar,f_cs,'r.-'); +# plt.plot(thetat,fact,'k-'); +# plt.title('Reconstruction with Compressed Sensing') #------------------------------------------------------------------------------------------------- @@ -396,7 +439,6 @@ # Calculate centered correlation f0 = factr - np.mean(factr); f1 = f_fourier - np.mean(f_fourier); - corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); f1 = f_capon - np.mean(f_capon); corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); @@ -499,11 +541,14 @@ ave = np.zeros(shape=(4)) print metric - -ave[0] = nanmean(metric[0,:,:]); # Fourier -ave[1] = nanmean(metric[1,:,:]); # Capon -ave[2] = nanmean(metric[2,:,:]); # MaxEnt -ave[3] = nanmean(metric[3,:,:]); # Compressed Sensing +print metric.shape + +ave[0] = nanmean(metric[0,0,:]); # Fourier +ave[1] = nanmean(metric[1,0,:]); # Capon +ave[2] = nanmean(metric[2,0,:]); # MaxEnt +ave[3] = nanmean(metric[3,0,:]); # Compressed Sensing + +print ave # Plot based on chosen metric plt.figure(6); Index: trunk/src/src/dualfilt1.py =================================================================== diff --git a/trunk/src/src/dualfilt1.py b/trunk/src/src/dualfilt1.py --- a/trunk/src/src/dualfilt1.py (revision 25) +++ b/trunk/src/src/dualfilt1.py (revision 26) @@ -26,7 +26,7 @@ # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY # http://taco.poly.edu/WaveletSoftware/ - # These cofficients are rounded to 8 decimal places. + # These coefficients are rounded to 8 decimal places. a1 = np.array([ [ 0.03516384000000, 0], Index: trunk/src/src/irls_dn2.py =================================================================== diff --git a/trunk/src/src/irls_dn2.py b/trunk/src/src/irls_dn2.py --- a/trunk/src/src/irls_dn2.py (revision 25) +++ b/trunk/src/src/irls_dn2.py (revision 26) @@ -37,7 +37,7 @@ fid = norm(np.dot(A,u)-b)**2; print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid; - # Refinement using fzero/ brentq + # Refinement using brentq (equiv to fzero in matlab) lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing def myfun(lambda1):