|
|
@@
-1,527
+1,572
|
|
|
1
|
1
|
#!/usr/bin/env python
No newline at end of file
|
|
|
2
|
2
|
No newline at end of file
|
|
|
3
|
3
|
#----------------------------------------------------------
No newline at end of file
|
|
|
4
|
4
|
# Original MATLAB code developed by Brian Harding
No newline at end of file
|
|
|
5
|
5
|
# Rewritten in Python by Yolian Amaro
No newline at end of file
|
|
|
6
|
6
|
# Python version 2.7
No newline at end of file
|
|
|
7
|
7
|
# May 15, 2014
No newline at end of file
|
|
|
8
|
8
|
# Jicamarca Radio Observatory
No newline at end of file
|
|
|
9
|
9
|
#----------------------------------------------------------
No newline at end of file
|
|
|
10
|
10
|
No newline at end of file
|
|
|
11
|
11
|
import time
No newline at end of file
|
|
|
12
|
12
|
import matplotlib.pyplot as plt
No newline at end of file
|
|
|
13
|
13
|
from scipy.optimize import root
No newline at end of file
|
|
|
|
14
|
from scipy.optimize import newton_krylov
No newline at end of file
|
|
|
14
|
15
|
from scipy.stats import nanmean
No newline at end of file
|
|
|
15
|
16
|
No newline at end of file
|
|
|
16
|
17
|
from y_hysell96 import *
No newline at end of file
|
|
|
17
|
18
|
from deb4_basis import *
No newline at end of file
|
|
|
18
|
19
|
from modelf import *
No newline at end of file
|
|
|
19
|
20
|
from irls_dn2 import *
No newline at end of file
|
|
|
20
|
21
|
No newline at end of file
|
|
|
21
|
22
|
No newline at end of file
|
|
|
22
|
23
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
23
|
24
|
# Set parameters
No newline at end of file
|
|
|
24
|
25
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
25
|
26
|
No newline at end of file
|
|
|
26
|
27
|
## Calculate Forward Model
No newline at end of file
|
|
|
27
|
28
|
lambda1 = 6.0
No newline at end of file
|
|
|
28
|
29
|
k = 2*np.pi/lambda1
No newline at end of file
|
|
|
29
|
30
|
No newline at end of file
|
|
|
30
|
31
|
## Magnetic Declination
No newline at end of file
|
|
|
31
|
32
|
dec = -1.24
No newline at end of file
|
|
|
32
|
33
|
No newline at end of file
|
|
|
33
|
34
|
## Loads Jicamarca antenna positions
No newline at end of file
|
|
|
34
|
35
|
antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False)
No newline at end of file
|
|
|
35
|
36
|
No newline at end of file
|
|
|
36
|
37
|
## rx and ry -- for plotting purposes
No newline at end of file
|
|
|
37
|
38
|
rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
No newline at end of file
|
|
|
38
|
39
|
ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
No newline at end of file
|
|
|
39
|
40
|
No newline at end of file
|
|
|
40
|
41
|
## Plot of antenna positions
No newline at end of file
|
|
|
41
|
42
|
plt.figure(1)
No newline at end of file
|
|
|
42
|
43
|
plt.plot(rx, ry, 'ro')
No newline at end of file
|
|
|
43
|
44
|
plt.title("Jicamarca Antenna Positions")
No newline at end of file
|
|
|
44
|
45
|
plt.xlabel("rx")
No newline at end of file
|
|
|
45
|
46
|
plt.ylabel("ry")
No newline at end of file
|
|
|
46
|
47
|
plt.grid(True)
No newline at end of file
|
|
|
47
|
48
|
plt.plot(rx, rx, 'b-')
No newline at end of file
|
|
|
48
|
49
|
plt.draw()
No newline at end of file
|
|
|
49
|
50
|
No newline at end of file
|
|
|
50
|
51
|
## Jicamarca is nominally at a 45 degree angle
No newline at end of file
|
|
|
51
|
52
|
theta = 45 - dec;
No newline at end of file
|
|
|
52
|
53
|
No newline at end of file
|
|
|
53
|
54
|
## Rotation matrix from antenna coord to magnetic coord (East North) -> 45 degrees
No newline at end of file
|
|
|
54
|
55
|
theta_rad = np.radians(theta) # trig functions take radians as argument
No newline at end of file
|
|
|
55
|
56
|
val1 = float( np.cos(theta_rad) )
|
|
|
|
57
|
No newline at end of file
|
|
|
56
|
|
val2 = float( np.sin(theta_rad) )
|
|
|
|
|
No newline at end of file
|
|
|
|
58
|
No newline at end of file
|
|
|
57
|
|
val3 = float( -1*np.sin(theta_rad))
No newline at end of file
|
|
|
58
|
59
|
val4 = float( np.cos(theta_rad) )
No newline at end of file
|
|
|
59
|
60
|
No newline at end of file
|
|
|
60
|
61
|
# Rotation matrix from antenna coord to magnetic coord (East North)
|
|
|
|
62
|
No newline at end of file
|
|
|
61
|
|
R = np.array( [[val1, val3], [val2, val4]] )
No newline at end of file
|
|
|
62
|
63
|
No newline at end of file
|
|
|
63
|
64
|
# Rotate antenna positions to magnetic coord.
No newline at end of file
|
|
|
64
|
65
|
AR = np.dot(R.T, antpos)
No newline at end of file
|
|
|
65
|
66
|
|
|
|
|
67
|
No newline at end of file
|
|
|
66
|
|
plt.plot(AR[0], 0*AR[1], 'yo')
No newline at end of file
|
|
|
67
|
68
|
plt.title("Jicamarca Antenna Positions")
No newline at end of file
|
|
|
68
|
69
|
plt.xlabel("rx")
No newline at end of file
|
|
|
69
|
70
|
plt.ylabel("ry")
No newline at end of file
|
|
|
70
|
71
|
plt.grid(True)
No newline at end of file
|
|
|
71
|
72
|
plt.plot(rx, 0*rx, 'g-')
No newline at end of file
|
|
|
72
|
73
|
plt.draw()
No newline at end of file
|
|
|
73
|
74
|
|
|
|
|
|
No newline at end of file
|
|
|
74
|
|
print antpos
|
|
|
|
|
No newline at end of file
|
|
|
75
|
|
print "\n"
|
|
|
|
|
No newline at end of file
|
|
|
76
|
|
print AR
No newline at end of file
|
|
|
77
|
75
|
No newline at end of file
|
|
|
78
|
76
|
# Only take the East component
No newline at end of file
|
|
|
79
|
77
|
r = AR[0,:]
No newline at end of file
|
|
|
80
|
78
|
r.sort()
|
|
|
|
|
No newline at end of file
|
|
|
81
|
|
|
|
|
|
|
No newline at end of file
|
|
|
82
|
|
# ---NEW---
|
|
|
|
|
No newline at end of file
|
|
|
83
|
|
r1 = AR[1,:] # longitudes / distances
|
|
|
|
|
No newline at end of file
|
|
|
84
|
|
d = 40 # separation of antennas in meters
|
|
|
|
|
No newline at end of file
|
|
|
85
|
|
theta_w = 15 # angle between antenna and point of observation in degrees
|
|
|
|
|
No newline at end of file
|
|
|
86
|
|
# arg = k*(r + d*np.cos(theta_w)) # this is the exponent argument in the corr integral
|
|
|
|
|
No newline at end of file
|
|
|
87
|
|
# --------
No newline at end of file
|
|
|
88
|
79
|
No newline at end of file
|
|
|
89
|
80
|
No newline at end of file
|
|
|
90
|
81
|
# Truth model (high and low resolution)
No newline at end of file
|
|
|
91
|
82
|
Nt = (1024.0)*(16.0) # number of pixels in truth image: high resolution
No newline at end of file
|
|
|
92
|
83
|
thbound = 9.0/180*np.pi # the width of the domain in angle space
No newline at end of file
|
|
|
93
|
84
|
thetat = np.linspace(-thbound, thbound,Nt) # image domain
No newline at end of file
|
|
|
94
|
85
|
thetat = thetat.T # transpose
No newline at end of file
|
|
|
95
|
86
|
Nr = (256.0) # number of pixels in reconstructed image: low res
No newline at end of file
|
|
|
96
|
87
|
thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
No newline at end of file
|
|
|
97
|
88
|
thetar = thetar.T # transpose
No newline at end of file
|
|
|
98
|
89
|
No newline at end of file
|
|
|
99
|
90
|
No newline at end of file
|
|
|
100
|
91
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
101
|
92
|
# Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b.
No newline at end of file
|
|
|
102
|
93
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
103
|
94
|
No newline at end of file
|
|
|
104
|
95
|
# Triple Gaussian
No newline at end of file
|
|
|
105
|
96
|
# a = np.array([3, 5, 2]);
No newline at end of file
|
|
|
106
|
97
|
# mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi, 7.0/180*np.pi]);
No newline at end of file
|
|
|
107
|
98
|
# sig = np.array([2.0/180*np.pi, 1.5/180*np.pi, 0.3/180*np.pi]);
No newline at end of file
|
|
|
108
|
99
|
# b = 0; # background
No newline at end of file
|
|
|
109
|
100
|
No newline at end of file
|
|
|
110
|
101
|
# Double Gaussian
No newline at end of file
|
|
|
111
|
102
|
# a = np.array([3, 5]);
No newline at end of file
|
|
|
112
|
103
|
# mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi]);
No newline at end of file
|
|
|
113
|
104
|
# sig = np.array([2.0/180*np.pi, 1.5/180*np.pi]);
No newline at end of file
|
|
|
114
|
105
|
# b = 0; # background
No newline at end of file
|
|
|
115
|
106
|
No newline at end of file
|
|
|
116
|
107
|
# Single Gaussian
No newline at end of file
|
|
|
117
|
108
|
a = np.array( [3] ) # amplitude(s)
No newline at end of file
|
|
|
118
|
109
|
mu = np.array( [-3.0/180*np.pi] ) # center
No newline at end of file
|
|
|
119
|
110
|
sig = np.array( [2.0/180*np.pi] ) # width
No newline at end of file
|
|
|
120
|
111
|
b = 0 # background constant
No newline at end of file
|
|
|
121
|
112
|
No newline at end of file
|
|
|
122
|
113
|
# Empty matrices for factors
No newline at end of file
|
|
|
123
|
114
|
fact = np.zeros(shape=(Nt,1))
No newline at end of file
|
|
|
124
|
115
|
factr = np.zeros(shape=(Nr,1))
No newline at end of file
|
|
|
125
|
116
|
|
|
|
|
117
|
No newline at end of file
|
|
|
126
|
|
# Constructing Gaussian model
No newline at end of file
|
|
|
127
|
118
|
for i in range(0, a.size):
No newline at end of file
|
|
|
128
|
119
|
temp = (-(thetat-mu[i])**2/(sig[i]**2))
No newline at end of file
|
|
|
129
|
120
|
tempr = (-(thetar-mu[i])**2/(sig[i]**2))
No newline at end of file
|
|
|
130
|
121
|
for j in range(0, temp.size):
No newline at end of file
|
|
|
131
|
122
|
fact[j] = fact[j] + a[i]*np.exp(temp[j]);
No newline at end of file
|
|
|
132
|
123
|
for m in range(0, tempr.size):
No newline at end of file
|
|
|
133
|
124
|
factr[m] = factr[m] + a[i]*np.exp(tempr[m]);
No newline at end of file
|
|
|
134
|
125
|
fact = fact + b;
No newline at end of file
|
|
|
135
|
126
|
factr = factr + b;
No newline at end of file
|
|
|
136
|
127
|
No newline at end of file
|
|
|
137
|
128
|
No newline at end of file
|
|
|
138
|
129
|
# #-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
139
|
130
|
# # Model for f: Square pulse
No newline at end of file
|
|
|
140
|
131
|
# #-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
141
|
132
|
# for j in range(0, fact.size):
No newline at end of file
|
|
|
142
|
133
|
# if (theta > -5.0/180*np.pi and theta < 2.0/180*np.pi):
No newline at end of file
|
|
|
143
|
134
|
# fact[j] = 0
No newline at end of file
|
|
|
144
|
135
|
# else:
No newline at end of file
|
|
|
145
|
136
|
# fact[j] = 1
No newline at end of file
|
|
|
146
|
137
|
# for k in range(0, factr.size):
No newline at end of file
|
|
|
147
|
138
|
# if (thetar[k] > -5.0/180*np.pi and thetar[k] < 2/180*np.pi):
No newline at end of file
|
|
|
148
|
139
|
# factr[k] = 0
No newline at end of file
|
|
|
149
|
140
|
# else:
No newline at end of file
|
|
|
150
|
141
|
# factr[k] = 1
No newline at end of file
|
|
|
151
|
142
|
No newline at end of file
|
|
|
152
|
143
|
# #-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
153
|
144
|
# # Model for f: Triangle pulse
No newline at end of file
|
|
|
154
|
145
|
# #-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
155
|
146
|
# mu = -1.0/180*np.pi;
No newline at end of file
|
|
|
156
|
147
|
# sig = 5.0/180*np.pi;
No newline at end of file
|
|
|
157
|
148
|
# wind1 = theta > mu-sig and theta < mu;
No newline at end of file
|
|
|
158
|
149
|
# wind2 = theta < mu+sig and theta > mu;
No newline at end of file
|
|
|
159
|
150
|
# fact = wind1 * (theta - (mu - sig));
No newline at end of file
|
|
|
160
|
151
|
# factr = wind1 * (thetar - (mu - sig));
No newline at end of file
|
|
|
161
|
152
|
# fact = fact + wind2 * (-(theta-(mu+sig)));
No newline at end of file
|
|
|
162
|
153
|
# factr = factr + wind2 * (-(thetar-(mu+sig)));
No newline at end of file
|
|
|
163
|
154
|
#
No newline at end of file
|
|
|
164
|
155
|
# fact = np.divide(fact, (np.sum(fact)*2*thbound/Nt)); # normalize to integral(f)==1
No newline at end of file
|
|
|
165
|
156
|
No newline at end of file
|
|
|
166
|
157
|
No newline at end of file
|
|
|
167
|
158
|
I = sum(fact)[0];
No newline at end of file
|
|
|
168
|
159
|
fact = fact/I; # normalize to sum(f)==1
No newline at end of file
|
|
|
169
|
160
|
factr = factr/I; # normalize to sum(f)==1
No newline at end of file
|
|
|
170
|
161
|
No newline at end of file
|
|
|
171
|
162
|
# Plot Gaussian pulse(s)
No newline at end of file
|
|
|
172
|
163
|
plt.figure(2)
|
|
|
|
164
|
No newline at end of file
|
|
|
173
|
|
plt.plot(thetat, fact, 'r--')
No newline at end of file
|
|
|
174
|
165
|
plt.plot(thetar, factr, 'ro')
|
|
|
|
166
|
No newline at end of file
|
|
|
175
|
|
plt.title("Truth models")
No newline at end of file
|
|
|
176
|
167
|
plt.ylabel("Amplitude")
|
|
|
|
168
|
No newline at end of file
|
|
|
177
|
|
plt.legend(('High res','Low res'), loc='upper right')
No newline at end of file
|
|
|
178
|
169
|
plt.draw()
No newline at end of file
|
|
|
179
|
170
|
No newline at end of file
|
|
|
180
|
171
|
No newline at end of file
|
|
|
181
|
172
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
182
|
173
|
# Control the type and number of inversions with:
No newline at end of file
|
|
|
183
|
174
|
# SNRdBvec: the SNRs that will be used.
No newline at end of file
|
|
|
184
|
175
|
# NN: the number of trials for each SNR
No newline at end of file
|
|
|
185
|
176
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
186
|
177
|
No newline at end of file
|
|
|
187
|
178
|
#SNRdBvec = np.linspace(5,20,10);
No newline at end of file
|
|
|
188
|
179
|
SNRdBvec = np.array([15]); # 15 dB
No newline at end of file
|
|
|
189
|
180
|
NN = 1; # number of trials at each SNR
No newline at end of file
|
|
|
190
|
181
|
No newline at end of file
|
|
|
191
|
182
|
# Statistics simulation (correlation, root mean square error)
No newline at end of file
|
|
|
192
|
183
|
corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
No newline at end of file
|
|
|
193
|
184
|
corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
No newline at end of file
|
|
|
194
|
185
|
rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
No newline at end of file
|
|
|
195
|
186
|
No newline at end of file
|
|
|
196
|
187
|
# For each SNR and trial
No newline at end of file
|
|
|
197
|
188
|
for snri in range(0, SNRdBvec.size):
No newline at end of file
|
|
|
198
|
189
|
SNRdB = SNRdBvec[snri];
No newline at end of file
|
|
|
199
|
190
|
SNR = 10**(SNRdB/10.0);
No newline at end of file
|
|
|
200
|
191
|
No newline at end of file
|
|
|
201
|
192
|
for Ni in range(0, NN):
No newline at end of file
|
|
|
|
193
|
|
|
|
|
|
No newline at end of file
|
|
|
|
194
|
# -------------------------------------- Voltages estimation----------------------------------------
|
|
|
|
|
No newline at end of file
|
|
|
|
195
|
v = np.zeros(shape=(r.size,1), dtype=object); # voltages matrix
|
|
|
|
|
No newline at end of file
|
|
|
|
196
|
R = np.zeros(shape=(r.size, r.size), dtype='complex128'); # corr matrix
|
|
|
|
|
No newline at end of file
|
|
|
|
197
|
mu1, sigma1 = 0, 1 # mean and standard deviation
|
|
|
|
|
No newline at end of file
|
|
|
|
198
|
num = 1
|
|
|
|
|
No newline at end of file
|
|
|
|
199
|
|
|
|
|
|
No newline at end of file
|
|
|
|
200
|
for t in range (0, 10):
|
|
|
|
|
No newline at end of file
|
|
|
|
201
|
num = t+1
|
|
|
|
|
No newline at end of file
|
|
|
|
202
|
print num
|
|
|
|
|
No newline at end of file
|
|
|
|
203
|
n1 = np.random.normal(mu1, sigma1, [Nt,1]) # random variable
|
|
|
|
|
No newline at end of file
|
|
|
|
204
|
n2 = np.random.normal(mu1, sigma1, [Nt,1]) # random variable
|
|
|
|
|
No newline at end of file
|
|
|
|
205
|
|
|
|
|
|
No newline at end of file
|
|
|
|
206
|
rho = (n1 + n2*1j)*np.sqrt(fact/2) # rho.size is Nt = 16834
|
|
|
|
|
No newline at end of file
|
|
|
|
207
|
# fact is size (Nt, 1) -> column vector
|
|
|
|
|
No newline at end of file
|
|
|
|
208
|
thetat = thetat.reshape(Nt,1) # thetat is now (Nt, 1)
|
|
|
|
|
No newline at end of file
|
|
|
|
209
|
|
|
|
|
|
No newline at end of file
|
|
|
|
210
|
for i in range(0, r.size):
|
|
|
|
|
No newline at end of file
|
|
|
|
211
|
v[i] = sum((rho)*np.exp(-1j*k*((r[i])*np.cos(thetat)))); # v is sum ( (Nt, 1)) = 1 value
|
|
|
|
|
No newline at end of file
|
|
|
|
212
|
|
|
|
|
|
No newline at end of file
|
|
|
|
213
|
# Computing correlations
|
|
|
|
|
No newline at end of file
|
|
|
|
214
|
for i1 in range(0, r.size):
|
|
|
|
|
No newline at end of file
|
|
|
|
215
|
for i2 in range(0,r.size):
|
|
|
|
|
No newline at end of file
|
|
|
|
216
|
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));
|
|
|
|
|
No newline at end of file
|
|
|
|
217
|
|
|
|
|
|
No newline at end of file
|
|
|
|
218
|
R = R / num
|
|
|
|
|
No newline at end of file
|
|
|
|
219
|
print "R from voltages"
|
|
|
|
|
No newline at end of file
|
|
|
|
220
|
print R
|
|
|
|
|
No newline at end of file
|
|
|
|
221
|
|
|
|
|
|
No newline at end of file
|
|
|
|
222
|
plt.figure(3);
|
|
|
|
|
No newline at end of file
|
|
|
|
223
|
plt.pcolor(abs(R));
|
|
|
|
|
No newline at end of file
|
|
|
|
224
|
plt.colorbar();
|
|
|
|
|
No newline at end of file
|
|
|
|
225
|
plt.title("Antenna Correlations")
|
|
|
|
|
No newline at end of file
|
|
|
|
226
|
plt.draw()
|
|
|
|
|
No newline at end of file
|
|
|
|
227
|
plt.show()
|
|
|
|
|
No newline at end of file
|
|
|
|
228
|
|
|
|
|
|
No newline at end of file
|
|
|
|
229
|
#--------------------------------------------------------------------------------------------------
|
|
|
|
|
No newline at end of file
|
|
|
|
230
|
No newline at end of file
|
|
|
202
|
231
|
# Calculate cross-correlation matrix (Fourier components of image)
No newline at end of file
|
|
|
203
|
232
|
# This is an inefficient way to do this.
No newline at end of file
|
|
|
204
|
233
|
No newline at end of file
|
|
|
|
234
|
# MATLAB CORR::: R(i1,i2) = fact.' * exp(j*k*(r(i1)-r(i2))*sin(thetat));
|
|
|
|
|
No newline at end of file
|
|
|
|
235
|
No newline at end of file
|
|
|
205
|
236
|
R = np.zeros(shape=(r.size, r.size), dtype=object);
No newline at end of file
|
|
|
206
|
237
|
No newline at end of file
|
|
|
|
238
|
thetat = thetat.reshape(Nt,)
No newline at end of file
|
|
|
207
|
239
|
for i1 in range(0, r.size):
No newline at end of file
|
|
|
208
|
240
|
for i2 in range(0,r.size):
|
|
|
|
241
|
No newline at end of file
|
|
|
209
|
|
R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat))))
|
|
|
|
|
No newline at end of file
|
|
|
|
242
|
No newline at end of file
|
|
|
210
|
|
R[i1,i2] = sum(R[i1,i2])
No newline at end of file
|
|
|
|
|
No newline at end of file
|
|
|
|
243
|
print "R w/o voltages\n", R
|
|
|
|
|
No newline at end of file
|
|
|
|
244
|
# plt.figure(3);
|
|
|
|
|
No newline at end of file
|
|
|
|
245
|
# plt.pcolor(abs(R));
|
|
|
|
|
No newline at end of file
|
|
|
|
246
|
# plt.colorbar();
|
|
|
|
|
No newline at end of file
|
|
|
|
247
|
# plt.title("Antenna Correlations")
|
|
|
|
|
No newline at end of file
|
|
|
|
248
|
# plt.xlabel("Antenna")
|
|
|
|
|
No newline at end of file
|
|
|
|
249
|
# plt.ylabel("Antenna")
|
|
|
|
|
No newline at end of file
|
|
|
|
250
|
# plt.draw()
|
|
|
|
|
No newline at end of file
|
|
|
|
251
|
# plt.show()
|
|
|
|
|
No newline at end of file
|
|
|
|
252
|
No newline at end of file
|
|
|
211
|
253
|
No newline at end of file
|
|
|
212
|
254
|
# Add uncertainty
No newline at end of file
|
|
|
213
|
255
|
# This is an ad-hoc way of adding "noise". It models some combination of
No newline at end of file
|
|
|
214
|
256
|
# receiver noise and finite integration times. We could use a more
No newline at end of file
|
|
|
215
|
257
|
# advanced model (like in Yu et al 2000) in the future.
No newline at end of file
|
|
|
216
|
258
|
No newline at end of file
|
|
|
217
|
259
|
# This is a way of adding noise while maintaining the
No newline at end of file
|
|
|
218
|
260
|
# positive-semi-definiteness of the matrix.
No newline at end of file
|
|
|
219
|
261
|
No newline at end of file
|
|
|
220
|
262
|
U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
No newline at end of file
|
|
|
221
|
263
|
No newline at end of file
|
|
|
222
|
264
|
sigma_noise = (np.linalg.norm(U,'fro')/SNR);
No newline at end of file
|
|
|
223
|
265
|
No newline at end of file
|
|
|
224
|
266
|
nz = np.multiply((sigma_noise * (np.random.randn(U.shape[0]) + 1j*np.random.randn(U.shape[0]))/np.sqrt(2)), (abs(U) > 0).astype(float));
No newline at end of file
|
|
|
225
|
267
|
No newline at end of file
|
|
|
226
|
268
|
Unz = U + nz;
No newline at end of file
|
|
|
227
|
269
|
Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
|
|
|
|
270
|
No newline at end of file
|
|
|
228
|
|
plt.figure(3);
|
|
|
|
|
No newline at end of file
|
|
|
|
271
|
No newline at end of file
|
|
|
229
|
|
plt.pcolor(abs(R));
|
|
|
|
|
No newline at end of file
|
|
|
|
272
|
No newline at end of file
|
|
|
230
|
|
plt.colorbar();
No newline at end of file
|
|
|
231
|
273
|
No newline at end of file
|
|
|
232
|
274
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
233
|
275
|
# Fourier Inversion
No newline at end of file
|
|
|
234
|
276
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
235
|
277
|
f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
No newline at end of file
|
|
|
236
|
278
|
No newline at end of file
|
|
|
237
|
279
|
for i in range(0, thetar.size):
No newline at end of file
|
|
|
238
|
280
|
th = thetar[i];
No newline at end of file
|
|
|
239
|
281
|
w = np.exp(1j*k*np.dot(r,np.sin(th)));
|
|
|
|
282
|
No newline at end of file
|
|
|
240
|
|
temp = np.dot(w.T.conj(),U)
|
|
|
|
|
No newline at end of file
|
|
|
241
|
|
f_fourier[i] = np.dot(temp, w);
No newline at end of file
|
|
|
242
|
283
|
No newline at end of file
|
|
|
243
|
284
|
f_fourier = f_fourier.real; # get rid of numerical imaginary noise
No newline at end of file
|
|
|
244
|
285
|
No newline at end of file
|
|
|
245
|
286
|
No newline at end of file
|
|
|
246
|
287
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
247
|
288
|
# Capon Inversion
No newline at end of file
|
|
|
248
|
289
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
249
|
290
|
f_capon = np.zeros(shape=(Nr,1));
No newline at end of file
|
|
|
250
|
291
|
No newline at end of file
|
|
|
251
|
292
|
tic_capon = time.time();
No newline at end of file
|
|
|
252
|
293
|
No newline at end of file
|
|
|
253
|
294
|
for i in range(0, thetar.size):
No newline at end of file
|
|
|
254
|
295
|
th = thetar[i];
No newline at end of file
|
|
|
255
|
296
|
w = np.exp(1j*k*np.dot(r,np.sin(th)));
No newline at end of file
|
|
|
256
|
297
|
f_capon[i] = 1/ ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real
No newline at end of file
|
|
|
257
|
298
|
No newline at end of file
|
|
|
258
|
299
|
toc_capon = time.time()
No newline at end of file
|
|
|
259
|
300
|
elapsed_time_capon = toc_capon - tic_capon;
No newline at end of file
|
|
|
260
|
301
|
No newline at end of file
|
|
|
261
|
302
|
f_capon = f_capon.real; # get rid of numerical imaginary noise
No newline at end of file
|
|
|
262
|
303
|
No newline at end of file
|
|
|
263
|
304
|
No newline at end of file
|
|
|
264
|
305
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
265
|
306
|
# MaxEnt Inversion
No newline at end of file
|
|
|
266
|
307
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
267
|
308
|
No newline at end of file
|
|
|
268
|
309
|
# Create the appropriate sensing matrix (split into real and imaginary # parts)
No newline at end of file
|
|
|
269
|
310
|
M = (r.size-1)*(r.size);
No newline at end of file
|
|
|
270
|
311
|
Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
No newline at end of file
|
|
|
271
|
312
|
Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
No newline at end of file
|
|
|
272
|
313
|
No newline at end of file
|
|
|
273
|
314
|
# Need to re-index our measurements from matrix R into vector g
No newline at end of file
|
|
|
274
|
315
|
g = np.zeros(shape=(M,1));
No newline at end of file
|
|
|
275
|
316
|
gnz = np.zeros(shape=(M,1)); # noisy version of g
No newline at end of file
|
|
|
276
|
317
|
No newline at end of file
|
|
|
277
|
318
|
# Triangular indexing to perform this re-indexing
No newline at end of file
|
|
|
278
|
319
|
T = np.ones(shape=(r.size,r.size));
No newline at end of file
|
|
|
279
|
320
|
[i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
No newline at end of file
|
|
|
280
|
321
|
No newline at end of file
|
|
|
281
|
322
|
# Build H
No newline at end of file
|
|
|
282
|
323
|
for i1 in range(0, r.size):
No newline at end of file
|
|
|
283
|
324
|
for i2 in range(i1+1, r.size):
|
|
|
|
325
|
No newline at end of file
|
|
|
284
|
|
idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward 1->27
No newline at end of file
|
|
|
285
|
326
|
idx1 = 2*idx;
No newline at end of file
|
|
|
286
|
327
|
idx2 = 2*idx+1;
No newline at end of file
|
|
|
287
|
328
|
Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj();
No newline at end of file
|
|
|
288
|
329
|
Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj();
No newline at end of file
|
|
|
289
|
330
|
Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt;
No newline at end of file
|
|
|
290
|
331
|
Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt;
No newline at end of file
|
|
|
291
|
332
|
g[idx1] = (R[i1,i2]).real*Nr/Nt;
No newline at end of file
|
|
|
292
|
333
|
g[idx2] = (R[i1,i2]).imag*Nr/Nt;
No newline at end of file
|
|
|
293
|
334
|
gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
No newline at end of file
|
|
|
294
|
335
|
gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
No newline at end of file
|
|
|
295
|
336
|
No newline at end of file
|
|
|
296
|
337
|
# Inversion
No newline at end of file
|
|
|
297
|
338
|
F = Nr/Nt; # normalization
No newline at end of file
|
|
|
298
|
339
|
sigma = 1; # set to 1 because the difference is accounted for in G
No newline at end of file
|
|
|
299
|
340
|
No newline at end of file
|
|
|
300
|
341
|
G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2
No newline at end of file
|
|
|
301
|
342
|
lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
No newline at end of file
|
|
|
302
|
343
|
No newline at end of file
|
|
|
303
|
344
|
No newline at end of file
|
|
|
304
|
345
|
# Whitened solution
No newline at end of file
|
|
|
305
|
346
|
def myfun(lambda1):
No newline at end of file
|
|
|
306
|
347
|
return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
No newline at end of file
|
|
|
307
|
348
|
No newline at end of file
|
|
|
308
|
349
|
No newline at end of file
|
|
|
309
|
350
|
tic_maxEnt = time.time(); # start time maxEnt
No newline at end of file
|
|
|
310
|
351
|
|
|
|
|
352
|
No newline at end of file
|
|
|
311
|
|
lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
No newline at end of file
|
|
|
|
|
No newline at end of file
|
|
|
|
353
|
# lambda1 = newton_krylov(myfun,lambda0, method='lgmres', x_tol=1e-14)
No newline at end of file
|
|
|
312
|
354
|
No newline at end of file
|
|
|
313
|
355
|
toc_maxEnt = time.time()
No newline at end of file
|
|
|
314
|
356
|
elapsed_time_maxent = toc_maxEnt - tic_maxEnt;
No newline at end of file
|
|
|
315
|
357
|
No newline at end of file
|
|
|
316
|
358
|
# Solution
No newline at end of file
|
|
|
317
|
359
|
lambda1 = lambda1.x;
No newline at end of file
|
|
|
318
|
360
|
No newline at end of file
|
|
|
319
|
361
|
f_maxent = modelf(lambda1, Hr, F);
No newline at end of file
|
|
|
320
|
362
|
No newline at end of file
|
|
|
321
|
363
|
#------ what is this needed for?-------
|
|
|
|
364
|
No newline at end of file
|
|
|
322
|
|
ystar = myfun(lambda1);
|
|
|
|
|
No newline at end of file
|
|
|
|
365
|
No newline at end of file
|
|
|
323
|
|
Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G));
|
|
|
|
|
No newline at end of file
|
|
|
|
366
|
No newline at end of file
|
|
|
324
|
|
ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
|
|
|
|
|
No newline at end of file
|
|
|
|
367
|
No newline at end of file
|
|
|
325
|
|
es = np.dot(Hr, f_maxent) - gnz;
|
|
|
|
|
No newline at end of file
|
|
|
|
368
|
No newline at end of file
|
|
|
326
|
|
chi2 = np.sum((es/sigma)**2);
No newline at end of file
|
|
|
327
|
369
|
#--------------------------------------
No newline at end of file
|
|
|
328
|
370
|
No newline at end of file
|
|
|
329
|
371
|
No newline at end of file
|
|
|
330
|
372
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
331
|
373
|
# CS inversion using Iteratively Reweighted Least Squares (IRLS)
No newline at end of file
|
|
|
332
|
374
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
333
|
375
|
# (Use Nr, thetar, gnz, and Hr from MaxEnt above)
No newline at end of file
|
|
|
334
|
376
|
No newline at end of file
|
|
|
335
|
377
|
Psi = deb4_basis(Nr);
No newline at end of file
|
|
|
336
|
378
|
No newline at end of file
|
|
|
337
|
379
|
# add "sum to 1" constraint
No newline at end of file
|
|
|
338
|
380
|
H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
No newline at end of file
|
|
|
339
|
381
|
g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
No newline at end of file
|
|
|
340
|
382
|
No newline at end of file
|
|
|
341
|
383
|
tic_cs = time.time();
No newline at end of file
|
|
|
342
|
384
|
No newline at end of file
|
|
|
343
|
385
|
# Inversion
No newline at end of file
|
|
|
344
|
386
|
s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
No newline at end of file
|
|
|
345
|
387
|
No newline at end of file
|
|
|
346
|
388
|
toc_cs = time.time()
No newline at end of file
|
|
|
347
|
389
|
elapsed_time_cs = toc_cs - tic_cs;
No newline at end of file
|
|
|
348
|
390
|
No newline at end of file
|
|
|
349
|
391
|
# Brightness function
No newline at end of file
|
|
|
350
|
392
|
f_cs = np.dot(Psi,s);
No newline at end of file
|
|
|
351
|
393
|
No newline at end of file
|
|
|
352
|
394
|
No newline at end of file
|
|
|
353
|
395
|
print "shapes:"
No newline at end of file
|
|
|
354
|
396
|
print "Psi ", Psi.shape
No newline at end of file
|
|
|
355
|
397
|
print "s ", s.shape
No newline at end of file
|
|
|
356
|
398
|
print "f_CS ", f_cs.shape
No newline at end of file
|
|
|
357
|
399
|
|
|
|
|
400
|
No newline at end of file
|
|
|
358
|
|
# Plot
|
|
|
|
|
No newline at end of file
|
|
|
|
401
|
No newline at end of file
|
|
|
359
|
|
plt.figure(4)
|
|
|
|
|
No newline at end of file
|
|
|
|
402
|
No newline at end of file
|
|
|
360
|
|
plt.plot(thetar,f_cs,'r.-');
|
|
|
|
|
No newline at end of file
|
|
|
|
403
|
No newline at end of file
|
|
|
361
|
|
plt.plot(thetat,fact,'k-');
|
|
|
|
|
No newline at end of file
|
|
|
|
404
|
No newline at end of file
|
|
|
362
|
|
plt.title('Reconstruction with Compressed Sensing')
No newline at end of file
|
|
|
|
|
No newline at end of file
|
|
|
|
405
|
# plt.title('Reconstruction with Compressed Sensing')
No newline at end of file
|
|
|
363
|
406
|
No newline at end of file
|
|
|
364
|
407
|
No newline at end of file
|
|
|
365
|
408
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
366
|
409
|
# Scaling and shifting
No newline at end of file
|
|
|
367
|
410
|
# (Only necessary for Capon solution)
No newline at end of file
|
|
|
368
|
411
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
369
|
412
|
f_capon = f_capon/np.max(f_capon)*np.max(fact);
No newline at end of file
|
|
|
370
|
413
|
No newline at end of file
|
|
|
371
|
414
|
No newline at end of file
|
|
|
372
|
415
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
373
|
416
|
# Analyze stuff
No newline at end of file
|
|
|
374
|
417
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
375
|
418
|
No newline at end of file
|
|
|
376
|
419
|
# Calculate MSE
No newline at end of file
|
|
|
377
|
420
|
rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
No newline at end of file
|
|
|
378
|
421
|
rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
No newline at end of file
|
|
|
379
|
422
|
rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
No newline at end of file
|
|
|
380
|
423
|
rmse_cs = np.sqrt(np.mean((f_cs - factr)**2));
No newline at end of file
|
|
|
381
|
424
|
No newline at end of file
|
|
|
382
|
425
|
No newline at end of file
|
|
|
383
|
426
|
relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
No newline at end of file
|
|
|
384
|
427
|
relrmse_capon = rmse_capon / np.linalg.norm(fact);
No newline at end of file
|
|
|
385
|
428
|
relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
No newline at end of file
|
|
|
386
|
429
|
relrmse_cs = rmse_cs / np.linalg.norm(fact);
No newline at end of file
|
|
|
387
|
430
|
No newline at end of file
|
|
|
388
|
431
|
No newline at end of file
|
|
|
389
|
432
|
# Calculate correlation
No newline at end of file
|
|
|
390
|
433
|
corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
No newline at end of file
|
|
|
391
|
434
|
corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
No newline at end of file
|
|
|
392
|
435
|
corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
No newline at end of file
|
|
|
393
|
436
|
corr_cs = np.dot(f_cs.T.conj(),factr) / (np.linalg.norm(f_cs)*np.linalg.norm(factr));
No newline at end of file
|
|
|
394
|
437
|
No newline at end of file
|
|
|
395
|
438
|
No newline at end of file
|
|
|
396
|
439
|
# Calculate centered correlation
No newline at end of file
|
|
|
397
|
440
|
f0 = factr - np.mean(factr);
No newline at end of file
|
|
|
398
|
441
|
f1 = f_fourier - np.mean(f_fourier);
No newline at end of file
|
|
|
399
|
442
|
No newline at end of file
|
|
|
400
|
443
|
corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
No newline at end of file
|
|
|
401
|
444
|
f1 = f_capon - np.mean(f_capon);
No newline at end of file
|
|
|
402
|
445
|
corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
No newline at end of file
|
|
|
403
|
446
|
f1 = f_maxent - np.mean(f_maxent);
No newline at end of file
|
|
|
404
|
447
|
corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
No newline at end of file
|
|
|
405
|
448
|
f1 = f_cs - np.mean(f_cs);
No newline at end of file
|
|
|
406
|
449
|
corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
No newline at end of file
|
|
|
407
|
450
|
No newline at end of file
|
|
|
408
|
451
|
No newline at end of file
|
|
|
409
|
452
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
410
|
453
|
# Plot stuff
No newline at end of file
|
|
|
411
|
454
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
412
|
455
|
No newline at end of file
|
|
|
413
|
456
|
#---- Capon----#
No newline at end of file
|
|
|
414
|
457
|
plt.figure(5)
No newline at end of file
|
|
|
415
|
458
|
plt.subplot(3, 1, 1)
No newline at end of file
|
|
|
416
|
459
|
plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon')
No newline at end of file
|
|
|
417
|
460
|
plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
No newline at end of file
|
|
|
418
|
461
|
plt.ylabel('Power (arbitrary units)')
No newline at end of file
|
|
|
419
|
462
|
plt.legend(loc='upper right')
No newline at end of file
|
|
|
420
|
463
|
No newline at end of file
|
|
|
421
|
464
|
# formatting y-axis
No newline at end of file
|
|
|
422
|
465
|
locs,labels = plt.yticks()
No newline at end of file
|
|
|
423
|
466
|
plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
No newline at end of file
|
|
|
424
|
467
|
plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
No newline at end of file
|
|
|
425
|
468
|
No newline at end of file
|
|
|
426
|
469
|
No newline at end of file
|
|
|
427
|
470
|
#---- MaxEnt----#
No newline at end of file
|
|
|
428
|
471
|
plt.subplot(3, 1, 2)
No newline at end of file
|
|
|
429
|
472
|
plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt')
No newline at end of file
|
|
|
430
|
473
|
plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
No newline at end of file
|
|
|
431
|
474
|
plt.ylabel('Power (arbitrary units)')
No newline at end of file
|
|
|
432
|
475
|
plt.legend(loc='upper right')
No newline at end of file
|
|
|
433
|
476
|
No newline at end of file
|
|
|
434
|
477
|
# formatting y-axis
No newline at end of file
|
|
|
435
|
478
|
locs,labels = plt.yticks()
No newline at end of file
|
|
|
436
|
479
|
plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
No newline at end of file
|
|
|
437
|
480
|
plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
No newline at end of file
|
|
|
438
|
481
|
No newline at end of file
|
|
|
439
|
482
|
No newline at end of file
|
|
|
440
|
483
|
#---- Compressed Sensing----#
No newline at end of file
|
|
|
441
|
484
|
plt.subplot(3, 1, 3)
No newline at end of file
|
|
|
442
|
485
|
plt.plot(180/np.pi*thetar, f_cs, 'r', label='CS')
No newline at end of file
|
|
|
443
|
486
|
plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
No newline at end of file
|
|
|
444
|
487
|
plt.ylabel('Power (arbitrary units)')
No newline at end of file
|
|
|
445
|
488
|
plt.legend(loc='upper right')
No newline at end of file
|
|
|
446
|
489
|
No newline at end of file
|
|
|
447
|
490
|
# formatting y-axis
No newline at end of file
|
|
|
448
|
491
|
locs,labels = plt.yticks()
No newline at end of file
|
|
|
449
|
492
|
plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
No newline at end of file
|
|
|
450
|
493
|
plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
No newline at end of file
|
|
|
451
|
494
|
No newline at end of file
|
|
|
452
|
495
|
# # Store Results
No newline at end of file
|
|
|
453
|
496
|
corr[0, snri, Ni] = corr_fourier;
No newline at end of file
|
|
|
454
|
497
|
corr[1, snri, Ni] = corr_capon;
No newline at end of file
|
|
|
455
|
498
|
corr[2, snri, Ni] = corr_maxent;
No newline at end of file
|
|
|
456
|
499
|
corr[3, snri, Ni] = corr_cs;
No newline at end of file
|
|
|
457
|
500
|
No newline at end of file
|
|
|
458
|
501
|
rmse[0,snri,Ni] = relrmse_fourier;
No newline at end of file
|
|
|
459
|
502
|
rmse[1,snri,Ni] = relrmse_capon;
No newline at end of file
|
|
|
460
|
503
|
rmse[2,snri,Ni] = relrmse_maxent;
No newline at end of file
|
|
|
461
|
504
|
rmse[3,snri,Ni] = relrmse_cs;
No newline at end of file
|
|
|
462
|
505
|
No newline at end of file
|
|
|
463
|
506
|
corrc[0,snri,Ni] = corrc_fourier;
No newline at end of file
|
|
|
464
|
507
|
corrc[1,snri,Ni] = corrc_capon;
No newline at end of file
|
|
|
465
|
508
|
corrc[2,snri,Ni] = corrc_maxent;
No newline at end of file
|
|
|
466
|
509
|
corrc[3,snri,Ni] = corrc_cs;
No newline at end of file
|
|
|
467
|
510
|
No newline at end of file
|
|
|
468
|
511
|
print "--------Time performace--------"
No newline at end of file
|
|
|
469
|
512
|
print 'Capon:\t', elapsed_time_capon, 'sec';
No newline at end of file
|
|
|
470
|
513
|
print 'Maxent:\t',elapsed_time_maxent, 'sec';
No newline at end of file
|
|
|
471
|
514
|
print 'CS:\t',elapsed_time_cs, 'sec\n';
No newline at end of file
|
|
|
472
|
515
|
No newline at end of file
|
|
|
473
|
516
|
print (NN*(snri)+Ni+1), '/', (SNRdBvec.size*NN), '\n';
No newline at end of file
|
|
|
474
|
517
|
No newline at end of file
|
|
|
475
|
518
|
No newline at end of file
|
|
|
476
|
519
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
477
|
520
|
# Analyze and plot statistics
No newline at end of file
|
|
|
478
|
521
|
#-------------------------------------------------------------------------------------------------
No newline at end of file
|
|
|
479
|
522
|
No newline at end of file
|
|
|
480
|
523
|
metric = corr; # set this to rmse, corr, or corrc
No newline at end of file
|
|
|
481
|
524
|
No newline at end of file
|
|
|
482
|
525
|
# Remove outliers (this part was experimental and wasn't used in the paper)
No newline at end of file
|
|
|
483
|
526
|
# nsig = 3;
No newline at end of file
|
|
|
484
|
527
|
# for i = 1:4
No newline at end of file
|
|
|
485
|
528
|
# for snri = 1:length(SNRdBvec)
No newline at end of file
|
|
|
486
|
529
|
# av = mean(met(i,snri,:));
No newline at end of file
|
|
|
487
|
530
|
# s = std(met(i,snri,:));
No newline at end of file
|
|
|
488
|
531
|
# idx = abs(met(i,snri,:) - av) > nsig*s;
No newline at end of file
|
|
|
489
|
532
|
# met(i,snri,idx) = nan;
No newline at end of file
|
|
|
490
|
533
|
# if sum(idx)>0
No newline at end of file
|
|
|
491
|
534
|
# fprintf('i=%i, snr=%i, %i/%i pts removed\n',...
No newline at end of file
|
|
|
492
|
535
|
# i,round(SNRdBvec(snri)),sum(idx),length(idx));
No newline at end of file
|
|
|
493
|
536
|
# end
No newline at end of file
|
|
|
494
|
537
|
# end
No newline at end of file
|
|
|
495
|
538
|
# end
No newline at end of file
|
|
|
496
|
539
|
No newline at end of file
|
|
|
497
|
540
|
No newline at end of file
|
|
|
498
|
541
|
# Avg ignoring NaNs
No newline at end of file
|
|
|
499
|
542
|
ave = np.zeros(shape=(4))
No newline at end of file
|
|
|
500
|
543
|
No newline at end of file
|
|
|
501
|
544
|
print metric
|
|
|
|
545
|
No newline at end of file
|
|
|
502
|
|
|
|
|
|
|
No newline at end of file
|
|
|
|
546
|
No newline at end of file
|
|
|
503
|
|
ave[0] = nanmean(metric[0,:,:]); # Fourier
|
|
|
|
|
No newline at end of file
|
|
|
|
547
|
No newline at end of file
|
|
|
504
|
|
ave[1] = nanmean(metric[1,:,:]); # Capon
|
|
|
|
|
No newline at end of file
|
|
|
|
548
|
No newline at end of file
|
|
|
505
|
|
ave[2] = nanmean(metric[2,:,:]); # MaxEnt
|
|
|
|
|
No newline at end of file
|
|
|
|
549
|
No newline at end of file
|
|
|
506
|
|
ave[3] = nanmean(metric[3,:,:]); # Compressed Sensing
No newline at end of file
|
|
|
|
|
No newline at end of file
|
|
|
|
550
|
|
|
|
|
|
No newline at end of file
|
|
|
|
551
|
print ave
|
|
|
|
|
No newline at end of file
|
|
|
|
552
|
No newline at end of file
|
|
|
507
|
553
|
No newline at end of file
|
|
|
508
|
554
|
# Plot based on chosen metric
No newline at end of file
|
|
|
509
|
555
|
plt.figure(6);
No newline at end of file
|
|
|
510
|
556
|
f = plt.scatter(SNRdBvec, ave[0], marker='+', color='b', s=60); # Fourier
No newline at end of file
|
|
|
511
|
557
|
c = plt.scatter(SNRdBvec, ave[1], marker='o', color= 'g', s=60); # Capon
No newline at end of file
|
|
|
512
|
558
|
me= plt.scatter(SNRdBvec, ave[2], marker='s', color= 'c', s=60); # MaxEnt
No newline at end of file
|
|
|
513
|
559
|
cs= plt.scatter(SNRdBvec, ave[3], marker='*', color='r', s=60); # Compressed Sensing
No newline at end of file
|
|
|
514
|
560
|
No newline at end of file
|
|
|
515
|
561
|
plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right')
No newline at end of file
|
|
|
516
|
562
|
plt.xlabel('SNR')
No newline at end of file
|
|
|
517
|
563
|
plt.ylabel('Correlation with Truth')
No newline at end of file
|
|
|
518
|
564
|
No newline at end of file
|
|
|
519
|
565
|
print "--------Correlations--------"
No newline at end of file
|
|
|
520
|
566
|
print "Fourier:", ave[0]
No newline at end of file
|
|
|
521
|
567
|
print "Capon:\t", ave[1]
No newline at end of file
|
|
|
522
|
568
|
print "MaxEnt:\t", ave[2]
No newline at end of file
|
|
|
523
|
569
|
print "CS:\t", ave[3]
No newline at end of file
|
|
|
524
|
570
|
No newline at end of file
|
|
|
525
|
571
|
plt.show()
No newline at end of file
|
|
|
526
|
572
|
No newline at end of file
|
|
|
527
|
573
|
|
