|
|
@@
-1,325
+1,430
|
|
|
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 math
No newline at end of file
|
|
|
12
|
12
|
#import cmath
No newline at end of file
|
|
|
13
|
13
|
#import scipy
No newline at end of file
|
|
|
14
|
14
|
#import matplotlib
No newline at end of file
|
|
|
15
|
15
|
import numpy as np
No newline at end of file
|
|
|
16
|
16
|
#from numpy import linalg
No newline at end of file
|
|
|
17
|
17
|
import matplotlib.pyplot as plt
No newline at end of file
|
|
|
18
|
18
|
from scipy import linalg
No newline at end of file
|
|
|
19
|
19
|
import time
No newline at end of file
|
|
|
20
|
20
|
from y_hysell96 import*
No newline at end of file
|
|
|
21
|
21
|
from deb4_basis import *
No newline at end of file
|
|
|
22
|
22
|
from modelf import *
No newline at end of file
|
|
|
23
|
23
|
from scipy.optimize import fsolve
No newline at end of file
|
|
|
24
|
24
|
from scipy.optimize import root
No newline at end of file
|
|
|
|
25
|
import pywt
No newline at end of file
|
|
|
25
|
26
|
No newline at end of file
|
|
|
26
|
27
|
No newline at end of file
|
|
|
27
|
28
|
## Calculate Forward Model
No newline at end of file
|
|
|
28
|
29
|
lambda1 = 6.0
No newline at end of file
|
|
|
29
|
30
|
k = 2*math.pi/lambda1
No newline at end of file
|
|
|
30
|
31
|
No newline at end of file
|
|
|
31
|
32
|
## Calculate Magnetic Declination
No newline at end of file
|
|
|
32
|
33
|
No newline at end of file
|
|
|
33
|
34
|
# [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this
No newline at end of file
|
|
|
34
|
35
|
No newline at end of file
|
|
|
35
|
36
|
# or calculate it with the above function
No newline at end of file
|
|
|
36
|
37
|
dec = -1.24
No newline at end of file
|
|
|
37
|
38
|
No newline at end of file
|
|
|
38
|
39
|
# loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
No newline at end of file
|
|
|
39
|
40
|
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
|
|
|
40
|
41
|
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
|
|
|
41
|
42
|
No newline at end of file
|
|
|
42
|
43
|
antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240],
No newline at end of file
|
|
|
43
|
44
|
[127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
No newline at end of file
|
|
|
44
|
45
|
No newline at end of file
|
|
|
45
|
46
|
plt.figure(1)
No newline at end of file
|
|
|
46
|
47
|
plt.plot(rx, ry, 'ro')
No newline at end of file
|
|
|
47
|
48
|
plt.draw()
No newline at end of file
|
|
|
48
|
49
|
No newline at end of file
|
|
|
49
|
50
|
# Jicamarca is nominally at a 45 degree angle
No newline at end of file
|
|
|
50
|
51
|
theta = 45 - dec;
No newline at end of file
|
|
|
51
|
52
|
No newline at end of file
|
|
|
52
|
53
|
# Rotation matrix from antenna coord to magnetic coord (East North)
No newline at end of file
|
|
|
53
|
54
|
theta_rad = math.radians(theta) # trig functions take radians as argument
No newline at end of file
|
|
|
54
|
55
|
val1 = float( math.cos(theta_rad) )
No newline at end of file
|
|
|
55
|
56
|
val2 = float( math.sin(theta_rad) )
No newline at end of file
|
|
|
56
|
57
|
val3 = float( -1*math.sin(theta_rad))
No newline at end of file
|
|
|
57
|
58
|
val4 = float( math.cos(theta_rad) )
No newline at end of file
|
|
|
58
|
59
|
No newline at end of file
|
|
|
59
|
60
|
# Rotation matrix from antenna coord to magnetic coord (East North)
No newline at end of file
|
|
|
60
|
61
|
R = np.array( [[val1, val3], [val2, val4]] );
No newline at end of file
|
|
|
61
|
62
|
No newline at end of file
|
|
|
62
|
63
|
# Rotate antenna positions to magnetic coord.
No newline at end of file
|
|
|
63
|
64
|
AR = np.dot(R.T, antpos);
No newline at end of file
|
|
|
64
|
65
|
No newline at end of file
|
|
|
65
|
66
|
# Only take the East component
No newline at end of file
|
|
|
66
|
67
|
r = AR[0,:]
No newline at end of file
|
|
|
67
|
68
|
r.sort() # ROW VECTOR?
No newline at end of file
|
|
|
68
|
69
|
No newline at end of file
|
|
|
69
|
70
|
# Truth model (high and low resolution)
No newline at end of file
|
|
|
70
|
71
|
Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
No newline at end of file
|
|
|
71
|
72
|
thbound = 9.0/180*math.pi; # the width of the domain in angle space
No newline at end of file
|
|
|
72
|
73
|
thetat = np.linspace(-thbound, thbound,Nt) # image domain
No newline at end of file
|
|
|
73
|
74
|
thetat = np.transpose(thetat) # transpose # FUNCIONA??????????????????????????????
No newline at end of file
|
|
|
74
|
75
|
Nr = (256.0); # number of pixels in reconstructed image: low res
No newline at end of file
|
|
|
75
|
76
|
thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
No newline at end of file
|
|
|
76
|
77
|
thetar = np.transpose(thetar) #transpose # FUNCIONA?????????????????????????????
No newline at end of file
|
|
|
77
|
78
|
No newline at end of file
|
|
|
78
|
79
|
# Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
No newline at end of file
|
|
|
79
|
80
|
# background constant b.
No newline at end of file
|
|
|
80
|
81
|
No newline at end of file
|
|
|
81
|
82
|
# Triple Gaussian
No newline at end of file
|
|
|
82
|
83
|
# a = np.array([3, 5, 2]);
No newline at end of file
|
|
|
83
|
84
|
# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
No newline at end of file
|
|
|
84
|
85
|
# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
No newline at end of file
|
|
|
85
|
86
|
# b = 0; # background
No newline at end of file
|
|
|
86
|
87
|
No newline at end of file
|
|
|
87
|
88
|
# Double Gaussian
No newline at end of file
|
|
|
88
|
89
|
# a = np.array([3, 5]);
No newline at end of file
|
|
|
89
|
90
|
# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
No newline at end of file
|
|
|
90
|
91
|
# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
No newline at end of file
|
|
|
91
|
92
|
# b = 0; # background
No newline at end of file
|
|
|
92
|
93
|
No newline at end of file
|
|
|
93
|
94
|
# Single Gaussian
No newline at end of file
|
|
|
94
|
95
|
a = np.array( [3] );
No newline at end of file
|
|
|
95
|
96
|
mu = np.array( [-3.0/180*math.pi] )
No newline at end of file
|
|
|
96
|
97
|
sig = np.array( [2.0/180*math.pi] )
No newline at end of file
|
|
|
97
|
98
|
b = 0;
No newline at end of file
|
|
|
98
|
99
|
No newline at end of file
|
|
|
99
|
100
|
fact = np.zeros(shape=(Nt,1));
No newline at end of file
|
|
|
100
|
101
|
factr = np.zeros(shape=(Nr,1));
No newline at end of file
|
|
|
101
|
102
|
No newline at end of file
|
|
|
102
|
103
|
for i in range(0, a.size):
No newline at end of file
|
|
|
103
|
104
|
temp = (-(thetat-mu[i])**2/(sig[i]**2))
No newline at end of file
|
|
|
104
|
105
|
tempr = (-(thetar-mu[i])**2/(sig[i]**2))
No newline at end of file
|
|
|
105
|
106
|
for j in range(0, temp.size):
No newline at end of file
|
|
|
106
|
107
|
fact[j] = fact[j] + a[i]*math.exp(temp[j]);
No newline at end of file
|
|
|
107
|
108
|
for m in range(0, tempr.size):
No newline at end of file
|
|
|
108
|
109
|
factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
No newline at end of file
|
|
|
109
|
110
|
fact = fact + b;
No newline at end of file
|
|
|
110
|
111
|
factr = factr + b;
No newline at end of file
|
|
|
111
|
112
|
No newline at end of file
|
|
|
112
|
113
|
# # model for f: Square pulse
No newline at end of file
|
|
|
113
|
114
|
# for j in range(0, fact.size):
No newline at end of file
|
|
|
114
|
115
|
# if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi):
No newline at end of file
|
|
|
115
|
116
|
# fact[j] = 0
No newline at end of file
|
|
|
116
|
117
|
# else:
No newline at end of file
|
|
|
117
|
118
|
# fact[j] = 1
No newline at end of file
|
|
|
118
|
119
|
# for k in range(0, factr.size):
No newline at end of file
|
|
|
119
|
120
|
# if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
No newline at end of file
|
|
|
120
|
121
|
# fact[k] = 0
No newline at end of file
|
|
|
121
|
122
|
# else:
No newline at end of file
|
|
|
122
|
123
|
# fact[k] = 1
No newline at end of file
|
|
|
123
|
124
|
#
No newline at end of file
|
|
|
124
|
125
|
#
No newline at end of file
|
|
|
125
|
126
|
# # model for f: triangle pulse
No newline at end of file
|
|
|
126
|
127
|
# mu = -1.0/180*math.pi;
No newline at end of file
|
|
|
127
|
128
|
# sig = 5.0/180*math.pi;
No newline at end of file
|
|
|
128
|
129
|
# wind1 = theta > mu-sig and theta < mu;
No newline at end of file
|
|
|
129
|
130
|
# wind2 = theta < mu+sig and theta > mu;
No newline at end of file
|
|
|
130
|
131
|
# fact = wind1 * (theta - (mu - sig));
No newline at end of file
|
|
|
131
|
132
|
# factr = wind1 * (thetar - (mu - sig));
No newline at end of file
|
|
|
132
|
133
|
# fact = fact + wind2 * (-(theta-(mu+sig)));
No newline at end of file
|
|
|
133
|
134
|
# factr = factr + wind2 * (-(thetar-(mu+sig)));
No newline at end of file
|
|
|
134
|
135
|
No newline at end of file
|
|
|
135
|
136
|
|
|
|
|
137
|
No newline at end of file
|
|
|
136
|
|
# fact = fact/(sum(fact)[0]*2*thbound/Nt); % normalize to integral(f)==1
No newline at end of file
|
|
|
137
|
138
|
I = sum(fact)[0];
No newline at end of file
|
|
|
138
|
139
|
fact = fact/I; # normalize to sum(f)==1
No newline at end of file
|
|
|
139
|
140
|
factr = factr/I; # normalize to sum(f)==1
No newline at end of file
|
|
|
140
|
141
|
#plot(thetat,fact,'r'); hold on;
No newline at end of file
|
|
|
141
|
142
|
#plot(thetar,factr,'k.'); hold off;
No newline at end of file
|
|
|
142
|
143
|
#xlim([min(thetat) max(thetat)]);
No newline at end of file
|
|
|
143
|
144
|
No newline at end of file
|
|
|
144
|
145
|
#x = np.linspace(thetat.min(), thetat.max) ????
No newline at end of file
|
|
|
145
|
146
|
#for i in range(0, thetat.size):
No newline at end of file
|
|
|
146
|
147
|
plt.figure(2)
No newline at end of file
|
|
|
147
|
148
|
plt.plot(thetat, fact, 'r--')
No newline at end of file
|
|
|
148
|
149
|
plt.plot(thetar, factr, 'ro')
No newline at end of file
|
|
|
149
|
150
|
plt.draw()
No newline at end of file
|
|
|
150
|
151
|
# xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
No newline at end of file
|
|
|
151
|
152
|
No newline at end of file
|
|
|
152
|
153
|
No newline at end of file
|
|
|
153
|
154
|
##
No newline at end of file
|
|
|
154
|
155
|
# Control the type and number of inversions with:
No newline at end of file
|
|
|
155
|
156
|
# SNRdBvec: the SNRs that will be used.
No newline at end of file
|
|
|
156
|
157
|
# NN: the number of trials for each SNR
No newline at end of file
|
|
|
157
|
158
|
No newline at end of file
|
|
|
158
|
159
|
#SNRdBvec = np.linspace(5,20,10);
No newline at end of file
|
|
|
159
|
160
|
SNRdBvec = np.array([15]);
No newline at end of file
|
|
|
160
|
161
|
NN = 1; # number of trial at each SNR
No newline at end of file
|
|
|
161
|
162
|
No newline at end of file
|
|
|
162
|
163
|
# if using vector arguments should be: (4,SNRdBvec.size,NN)
No newline at end of file
|
|
|
163
|
164
|
corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
No newline at end of file
|
|
|
164
|
165
|
corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
No newline at end of file
|
|
|
165
|
166
|
rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
No newline at end of file
|
|
|
166
|
167
|
No newline at end of file
|
|
|
167
|
168
|
for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
No newline at end of file
|
|
|
168
|
169
|
for Ni in range(0, NN):
No newline at end of file
|
|
|
169
|
170
|
SNRdB = SNRdBvec[snri];
No newline at end of file
|
|
|
170
|
171
|
SNR = 10**(SNRdB/10.0);
No newline at end of file
|
|
|
171
|
172
|
No newline at end of file
|
|
|
172
|
173
|
# Calculate cross-correlation matrix (Fourier components of image)
No newline at end of file
|
|
|
173
|
174
|
# This is an inefficient way to do this.
No newline at end of file
|
|
|
174
|
175
|
R = np.zeros(shape=(r.size, r.size), dtype=object);
No newline at end of file
|
|
|
175
|
176
|
No newline at end of file
|
|
|
176
|
177
|
for i1 in range(0, r.size):
No newline at end of file
|
|
|
177
|
178
|
for i2 in range(0,r.size):
No newline at end of file
|
|
|
178
|
179
|
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
|
|
|
179
|
180
|
R[i1,i2] = sum(R[i1,i2])
No newline at end of file
|
|
|
180
|
181
|
No newline at end of file
|
|
|
181
|
182
|
# Add uncertainty
No newline at end of file
|
|
|
182
|
183
|
# This is an ad-hoc way of adding "noise". It models some combination of
No newline at end of file
|
|
|
183
|
184
|
# receiver noise and finite integration times. We could use a more
No newline at end of file
|
|
|
184
|
185
|
# advanced model (like in Yu et al 2000) in the future.
No newline at end of file
|
|
|
185
|
186
|
No newline at end of file
|
|
|
186
|
187
|
# This is a way of adding noise while maintaining the
No newline at end of file
|
|
|
187
|
188
|
# positive-semi-definiteness of the matrix.
No newline at end of file
|
|
|
188
|
189
|
No newline at end of file
|
|
|
189
|
190
|
U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
No newline at end of file
|
|
|
190
|
191
|
No newline at end of file
|
|
|
191
|
192
|
sigma_noise = (np.linalg.norm(U,'fro')/SNR);
No newline at end of file
|
|
|
192
|
193
|
No newline at end of file
|
|
|
193
|
194
|
temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
No newline at end of file
|
|
|
194
|
195
|
temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
No newline at end of file
|
|
|
195
|
196
|
temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
No newline at end of file
|
|
|
196
|
197
|
temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
No newline at end of file
|
|
|
197
|
198
|
No newline at end of file
|
|
|
198
|
199
|
nz = np.multiply(temp4, temp3)
No newline at end of file
|
|
|
199
|
200
|
No newline at end of file
|
|
|
200
|
201
|
#nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
No newline at end of file
|
|
|
201
|
202
|
No newline at end of file
|
|
|
202
|
203
|
No newline at end of file
|
|
|
203
|
204
|
#nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
No newline at end of file
|
|
|
204
|
205
|
#nz = np.dot(sigma_noise, (np.dot((np.random.rand(8,8) + j*np.random.rand(8,8))/math.sqrt(2.0) , (abs(U) > 0).astype(int))));
No newline at end of file
|
|
|
205
|
206
|
No newline at end of file
|
|
|
206
|
207
|
Unz = U + nz;
No newline at end of file
|
|
|
207
|
208
|
Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
No newline at end of file
|
|
|
208
|
209
|
plt.figure(3);
No newline at end of file
|
|
|
209
|
210
|
plt.pcolor(abs(Rnz));
No newline at end of file
|
|
|
210
|
211
|
plt.colorbar();
No newline at end of file
|
|
|
211
|
212
|
|
|
|
|
213
|
No newline at end of file
|
|
|
212
|
|
# Fourier Inversion %%%%%%%%%%%%%%%%%%%
No newline at end of file
|
|
|
213
|
214
|
f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
No newline at end of file
|
|
|
214
|
215
|
No newline at end of file
|
|
|
215
|
216
|
for i in range(0, thetar.size):
No newline at end of file
|
|
|
216
|
217
|
th = thetar[i];
No newline at end of file
|
|
|
217
|
218
|
w = np.exp(1j*k*np.dot(r,np.sin(th)));
No newline at end of file
|
|
|
218
|
219
|
No newline at end of file
|
|
|
219
|
220
|
temp = np.dot(w.T.conj(),U)
No newline at end of file
|
|
|
220
|
221
|
No newline at end of file
|
|
|
221
|
222
|
f_fourier[i] = np.dot(temp, w);
No newline at end of file
|
|
|
222
|
223
|
No newline at end of file
|
|
|
223
|
224
|
f_fourier = f_fourier.real; # get rid of numerical imaginary noise
No newline at end of file
|
|
|
224
|
225
|
No newline at end of file
|
|
|
225
|
226
|
#print f_fourier
No newline at end of file
|
|
|
226
|
227
|
No newline at end of file
|
|
|
227
|
228
|
|
|
|
|
229
|
No newline at end of file
|
|
|
228
|
|
# Capon Inversion %%%%%%%%%%%%%%%%%%%%%%
No newline at end of file
|
|
|
229
|
230
|
No newline at end of file
|
|
|
230
|
231
|
f_capon = np.zeros(shape=(Nr,1));
No newline at end of file
|
|
|
231
|
232
|
No newline at end of file
|
|
|
232
|
233
|
#tic_capon = time.time();
No newline at end of file
|
|
|
233
|
234
|
No newline at end of file
|
|
|
234
|
235
|
for i in range(0, thetar.size):
No newline at end of file
|
|
|
235
|
236
|
th = thetar[i];
No newline at end of file
|
|
|
236
|
237
|
w = np.exp(1j*k*np.dot(r,np.sin(th)));
No newline at end of file
|
|
|
237
|
238
|
f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
No newline at end of file
|
|
|
238
|
239
|
No newline at end of file
|
|
|
239
|
240
|
#toc_capon = time.time()
No newline at end of file
|
|
|
240
|
241
|
No newline at end of file
|
|
|
241
|
242
|
# elapsed_time_capon = toc_capon - tic_capon;
No newline at end of file
|
|
|
242
|
243
|
No newline at end of file
|
|
|
243
|
244
|
f_capon = f_capon.real; # get rid of numerical imaginary noise
No newline at end of file
|
|
|
244
|
245
|
|
|
|
|
246
|
No newline at end of file
|
|
|
245
|
|
# MaxEnt Inversion %%%%%%%%%%%%%%%%%%%%%
|
|
|
|
|
No newline at end of file
|
|
|
|
247
|
No newline at end of file
|
|
|
246
|
|
|
|
|
|
|
No newline at end of file
|
|
|
|
248
|
No newline at end of file
|
|
|
247
|
|
# create the appropriate sensing matrix (split into real and imaginary % parts)
No newline at end of file
|
|
|
248
|
249
|
M = (r.size-1)*(r.size);
No newline at end of file
|
|
|
249
|
250
|
Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
No newline at end of file
|
|
|
250
|
251
|
Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
No newline at end of file
|
|
|
251
|
252
|
No newline at end of file
|
|
|
252
|
253
|
# need to re-index our measurements from matrix R into vector g
No newline at end of file
|
|
|
253
|
254
|
g = np.zeros(shape=(M,1));
No newline at end of file
|
|
|
254
|
255
|
gnz = np.zeros(shape=(M,1)); # noisy version of g
No newline at end of file
|
|
|
255
|
256
|
No newline at end of file
|
|
|
256
|
257
|
# triangular indexing to perform this re-indexing
No newline at end of file
|
|
|
257
|
258
|
T = np.ones(shape=(r.size,r.size));
No newline at end of file
|
|
|
258
|
259
|
[i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
No newline at end of file
|
|
|
259
|
260
|
No newline at end of file
|
|
|
260
|
261
|
# build H
No newline at end of file
|
|
|
261
|
262
|
for i1 in range(0, r.size):
No newline at end of file
|
|
|
262
|
263
|
for i2 in range(i1+1, r.size):
No newline at end of file
|
|
|
263
|
264
|
idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward
No newline at end of file
|
|
|
264
|
265
|
idx1 = 2*idx; # because index starts at 0
No newline at end of file
|
|
|
265
|
266
|
idx2 = 2*idx+1;
No newline at end of file
|
|
|
266
|
267
|
Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T;
No newline at end of file
|
|
|
267
|
268
|
Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T;
No newline at end of file
|
|
|
268
|
269
|
Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
No newline at end of file
|
|
|
269
|
270
|
Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
No newline at end of file
|
|
|
270
|
271
|
g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later
No newline at end of file
|
|
|
271
|
272
|
g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again
No newline at end of file
|
|
|
272
|
273
|
gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
No newline at end of file
|
|
|
273
|
274
|
gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
No newline at end of file
|
|
|
274
|
275
|
No newline at end of file
|
|
|
275
|
276
|
# inversion
No newline at end of file
|
|
|
276
|
277
|
F = Nr/Nt; # normalization
No newline at end of file
|
|
|
277
|
278
|
sigma = 1; # set to 1 because the difference is accounted for in G
No newline at end of file
|
|
|
278
|
279
|
No newline at end of file
|
|
|
279
|
280
|
##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
No newline at end of file
|
|
|
280
|
281
|
G = np.linalg.norm(g-gnz)**2; # pretend we know in advance the actual value of chi^2
No newline at end of file
|
|
|
281
|
282
|
No newline at end of file
|
|
|
282
|
283
|
lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
No newline at end of file
|
|
|
283
|
284
|
No newline at end of file
|
|
|
284
|
285
|
# Whitened solution
No newline at end of file
|
|
|
285
|
286
|
def myfun(lambda1):
No newline at end of file
|
|
|
286
|
287
|
return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
No newline at end of file
|
|
|
287
|
288
|
No newline at end of file
|
|
|
288
|
289
|
tic_maxEnt = time.time();
No newline at end of file
|
|
|
289
|
290
|
No newline at end of file
|
|
|
290
|
291
|
#sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
No newline at end of file
|
|
|
291
|
292
|
lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
No newline at end of file
|
|
|
292
|
293
|
No newline at end of file
|
|
|
293
|
294
|
#print lambda1
No newline at end of file
|
|
|
294
|
295
|
#print lambda1.x
No newline at end of file
|
|
|
295
|
296
|
No newline at end of file
|
|
|
296
|
297
|
lambda1 = lambda1.x;
No newline at end of file
|
|
|
297
|
298
|
No newline at end of file
|
|
|
298
|
299
|
toc_maxEnt = time.time();
No newline at end of file
|
|
|
299
|
300
|
f_maxent = modelf(lambda1, Hr, F);
No newline at end of file
|
|
|
300
|
301
|
ystar = myfun(lambda1);
No newline at end of file
|
|
|
301
|
302
|
Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G));
No newline at end of file
|
|
|
302
|
303
|
ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
No newline at end of file
|
|
|
303
|
304
|
es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
No newline at end of file
|
|
|
304
|
305
|
chi2 = np.sum((es/sigma)**2);
No newline at end of file
|
|
|
305
|
306
|
No newline at end of file
|
|
|
306
|
307
|
No newline at end of file
|
|
|
307
|
308
|
|
|
|
|
309
|
No newline at end of file
|
|
|
308
|
|
# CS inversion using irls %%%%%%%%%%%%%%%%%%%%%%%%
No newline at end of file
|
|
|
309
|
310
|
No newline at end of file
|
|
|
310
|
311
|
# (Use Nr, thetar, gnz, and Hr from MaxEnt above)
No newline at end of file
|
|
|
311
|
312
|
|
|
|
|
313
|
No newline at end of file
|
|
|
312
|
|
Psi = deb4_basis(Nr);
|
|
|
|
|
No newline at end of file
|
|
|
|
314
|
No newline at end of file
|
|
|
313
|
|
|
|
|
|
|
No newline at end of file
|
|
|
|
315
|
No newline at end of file
|
|
|
314
|
|
# # add "sum to 1" constraint
|
|
|
|
|
No newline at end of file
|
|
|
|
316
|
No newline at end of file
|
|
|
315
|
|
# H2 = [Hr; ones(1,Nr)];
|
|
|
|
|
No newline at end of file
|
|
|
|
317
|
No newline at end of file
|
|
|
316
|
|
# g2 = [gnz; Nr/Nt];
|
|
|
|
|
No newline at end of file
|
|
|
|
318
|
No newline at end of file
|
|
|
317
|
|
#
|
|
|
|
|
No newline at end of file
|
|
|
|
319
|
No newline at end of file
|
|
|
318
|
|
# s = irls_dn2(H2*Psi,g2,0.5,G);
No newline at end of file
|
|
|
|
|
No newline at end of file
|
|
|
|
320
|
N_temp = np.array([[Nr/Nt]]);
|
|
|
|
|
No newline at end of file
|
|
|
|
321
|
g2 = np.concatenate( (gnz, N_temp), axis=0 );
|
|
|
|
|
No newline at end of file
|
|
|
|
322
|
|
|
|
|
|
No newline at end of file
|
|
|
|
323
|
|
|
|
|
|
No newline at end of file
|
|
|
|
324
|
#s = irls_dn2(H2*Psi,g2,0.5,G);
No newline at end of file
|
|
|
319
|
325
|
# f_cs = Psi*s;
No newline at end of file
|
|
|
320
|
326
|
#
No newline at end of file
|
|
|
321
|
327
|
# # plot
No newline at end of file
|
|
|
322
|
328
|
# plot(thetar,f_cs,'r.-');
No newline at end of file
|
|
|
323
|
329
|
# hold on;
No newline at end of file
|
|
|
324
|
330
|
# plot(thetat,fact,'k-');
No newline at end of file
|
|
|
325
|
331
|
# hold off;
No newline at end of file
|
|
|
|
332
|
|
|
|
|
|
No newline at end of file
|
|
|
|
333
|
|
|
|
|
|
No newline at end of file
|
|
|
|
334
|
# # # Scaling and shifting
|
|
|
|
|
No newline at end of file
|
|
|
|
335
|
# # # Only necessary for capon solution
|
|
|
|
|
No newline at end of file
|
|
|
|
336
|
|
|
|
|
|
No newline at end of file
|
|
|
|
337
|
|
|
|
|
|
No newline at end of file
|
|
|
|
338
|
f_capon = f_capon/np.max(f_capon)*np.max(fact);
|
|
|
|
|
No newline at end of file
|
|
|
|
339
|
|
|
|
|
|
No newline at end of file
|
|
|
|
340
|
|
|
|
|
|
No newline at end of file
|
|
|
|
341
|
### analyze stuff ######################
|
|
|
|
|
No newline at end of file
|
|
|
|
342
|
# calculate MSE
|
|
|
|
|
No newline at end of file
|
|
|
|
343
|
rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
|
|
|
|
|
No newline at end of file
|
|
|
|
344
|
rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
|
|
|
|
|
No newline at end of file
|
|
|
|
345
|
rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
|
|
|
|
|
No newline at end of file
|
|
|
|
346
|
#rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
|
|
|
|
|
No newline at end of file
|
|
|
|
347
|
|
|
|
|
|
No newline at end of file
|
|
|
|
348
|
relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
|
|
|
|
|
No newline at end of file
|
|
|
|
349
|
relrmse_capon = rmse_capon / np.linalg.norm(fact);
|
|
|
|
|
No newline at end of file
|
|
|
|
350
|
relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
|
|
|
|
|
No newline at end of file
|
|
|
|
351
|
#relrmse_cs = rmse_cs / np.norm(fact);
|
|
|
|
|
No newline at end of file
|
|
|
|
352
|
|
|
|
|
|
No newline at end of file
|
|
|
|
353
|
factr = factr.T.conj()
|
|
|
|
|
No newline at end of file
|
|
|
|
354
|
|
|
|
|
|
No newline at end of file
|
|
|
|
355
|
# calculate correlation
|
|
|
|
|
No newline at end of file
|
|
|
|
356
|
corr_fourier = np.dot(f_fourier,factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
|
|
|
|
|
No newline at end of file
|
|
|
|
357
|
corr_capon = np.dot(f_capon,factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
|
|
|
|
|
No newline at end of file
|
|
|
|
358
|
corr_maxent = np.dot(f_maxent,factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
|
|
|
|
|
No newline at end of file
|
|
|
|
359
|
#corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
|
|
|
|
|
No newline at end of file
|
|
|
|
360
|
|
|
|
|
|
No newline at end of file
|
|
|
|
361
|
# calculate centered correlation
|
|
|
|
|
No newline at end of file
|
|
|
|
362
|
f0 = factr - np.mean(factr);
|
|
|
|
|
No newline at end of file
|
|
|
|
363
|
f1 = f_fourier - np.mean(f_fourier);
|
|
|
|
|
No newline at end of file
|
|
|
|
364
|
corrc_fourier = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
|
|
|
|
|
No newline at end of file
|
|
|
|
365
|
f1 = f_capon - np.mean(f_capon);
|
|
|
|
|
No newline at end of file
|
|
|
|
366
|
corrc_capon = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
|
|
|
|
|
No newline at end of file
|
|
|
|
367
|
f1 = f_maxent - np.mean(f_maxent);
|
|
|
|
|
No newline at end of file
|
|
|
|
368
|
corrc_maxent = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
|
|
|
|
|
No newline at end of file
|
|
|
|
369
|
#f1 = f_cs - mean(f_cs);
|
|
|
|
|
No newline at end of file
|
|
|
|
370
|
#corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
|
|
|
|
|
No newline at end of file
|
|
|
|
371
|
|
|
|
|
|
No newline at end of file
|
|
|
|
372
|
|
|
|
|
|
No newline at end of file
|
|
|
|
373
|
|
|
|
|
|
No newline at end of file
|
|
|
|
374
|
# # # plot stuff #########################
|
|
|
|
|
No newline at end of file
|
|
|
|
375
|
|
|
|
|
|
No newline at end of file
|
|
|
|
376
|
#---- Capon----
|
|
|
|
|
No newline at end of file
|
|
|
|
377
|
plt.figure(4)
|
|
|
|
|
No newline at end of file
|
|
|
|
378
|
plt.subplot(2, 1, 1)
|
|
|
|
|
No newline at end of file
|
|
|
|
379
|
plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
|
|
|
|
|
No newline at end of file
|
|
|
|
380
|
plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
|
|
|
|
|
No newline at end of file
|
|
|
|
381
|
plt.ylabel('Power (arbitrary units)')
|
|
|
|
|
No newline at end of file
|
|
|
|
382
|
plt.legend(loc='upper right')
|
|
|
|
|
No newline at end of file
|
|
|
|
383
|
|
|
|
|
|
No newline at end of file
|
|
|
|
384
|
# formatting y-axis
|
|
|
|
|
No newline at end of file
|
|
|
|
385
|
locs,labels = plt.yticks()
|
|
|
|
|
No newline at end of file
|
|
|
|
386
|
plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
|
|
|
|
|
No newline at end of file
|
|
|
|
387
|
plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
|
|
|
|
|
No newline at end of file
|
|
|
|
388
|
|
|
|
|
|
No newline at end of file
|
|
|
|
389
|
|
|
|
|
|
No newline at end of file
|
|
|
|
390
|
#---- MaxEnt----
|
|
|
|
|
No newline at end of file
|
|
|
|
391
|
plt.subplot(2, 1, 2)
|
|
|
|
|
No newline at end of file
|
|
|
|
392
|
plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
|
|
|
|
|
No newline at end of file
|
|
|
|
393
|
plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
|
|
|
|
|
No newline at end of file
|
|
|
|
394
|
plt.ylabel('Power (arbitrary units)')
|
|
|
|
|
No newline at end of file
|
|
|
|
395
|
plt.legend(loc='upper right')
|
|
|
|
|
No newline at end of file
|
|
|
|
396
|
|
|
|
|
|
No newline at end of file
|
|
|
|
397
|
# formatting y-axis
|
|
|
|
|
No newline at end of file
|
|
|
|
398
|
locs,labels = plt.yticks()
|
|
|
|
|
No newline at end of file
|
|
|
|
399
|
plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
|
|
|
|
|
No newline at end of file
|
|
|
|
400
|
plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
|
|
|
|
|
No newline at end of file
|
|
|
|
401
|
|
|
|
|
|
No newline at end of file
|
|
|
|
402
|
plt.show()
|
|
|
|
|
No newline at end of file
|
|
|
|
403
|
|
|
|
|
|
No newline at end of file
|
|
|
|
404
|
|
|
|
|
|
No newline at end of file
|
|
|
|
405
|
# # subplot(3,1,2);
|
|
|
|
|
No newline at end of file
|
|
|
|
406
|
# # plot(180/pi*thetar,f_maxent,'r-');
|
|
|
|
|
No newline at end of file
|
|
|
|
407
|
# # hold on;
|
|
|
|
|
No newline at end of file
|
|
|
|
408
|
# # plot(180/pi*thetat,fact,'k--');
|
|
|
|
|
No newline at end of file
|
|
|
|
409
|
# # hold off;
|
|
|
|
|
No newline at end of file
|
|
|
|
410
|
# # ylim([min(f_cs) 1.1*max(fact)]);
|
|
|
|
|
No newline at end of file
|
|
|
|
411
|
# # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_maxent, corr_maxent, corrc_maxent));
|
|
|
|
|
No newline at end of file
|
|
|
|
412
|
# # ylabel({'Power';'(arbitrary units)'})
|
|
|
|
|
No newline at end of file
|
|
|
|
413
|
# # # title 'Maximum Entropy Method'
|
|
|
|
|
No newline at end of file
|
|
|
|
414
|
# # legend('MaxEnt','Truth');
|
|
|
|
|
No newline at end of file
|
|
|
|
415
|
# #
|
|
|
|
|
No newline at end of file
|
|
|
|
416
|
# # subplot(3,1,3);
|
|
|
|
|
No newline at end of file
|
|
|
|
417
|
# # plot(180/pi*thetar,f_cs,'r-');
|
|
|
|
|
No newline at end of file
|
|
|
|
418
|
# # hold on;
|
|
|
|
|
No newline at end of file
|
|
|
|
419
|
# # plot(180/pi*thetat,fact,'k--');
|
|
|
|
|
No newline at end of file
|
|
|
|
420
|
# # hold off;
|
|
|
|
|
No newline at end of file
|
|
|
|
421
|
# # ylim([min(f_cs) 1.1*max(fact)]);
|
|
|
|
|
No newline at end of file
|
|
|
|
422
|
# # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
|
|
|
|
|
No newline at end of file
|
|
|
|
423
|
# # # title 'Compressed Sensing - Debauchies Wavelets'
|
|
|
|
|
No newline at end of file
|
|
|
|
424
|
# # xlabel 'Degrees'
|
|
|
|
|
No newline at end of file
|
|
|
|
425
|
# # ylabel({'Power';'(arbitrary units)'})
|
|
|
|
|
No newline at end of file
|
|
|
|
426
|
# # legend('Comp. Sens.','Truth');
|
|
|
|
|
No newline at end of file
|
|
|
|
427
|
# #
|
|
|
|
|
No newline at end of file
|
|
|
|
428
|
# # # set(gcf,'Position',[749 143 528 881]); # CSL
|
|
|
|
|
No newline at end of file
|
|
|
|
429
|
# # # set(gcf,'Position',[885 -21 528 673]); # macbook
|
|
|
|
|
No newline at end of file
|
|
|
|
430
|
# # pause(0.01);
No newline at end of file
|
