##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
yamaro -
r21:22
parent child
Show More
Show file before | Show file after | Hide comments
@@ -2,127 +2,133
2 No newline at end of file
2
3 #---------------------------------------------------------- No newline at end of file
3 #----------------------------------------------------------
4 # Original MATLAB code developed by Brian Harding
4 # Original MATLAB code developed by Brian Harding
No newline at end of file
5 # Rewritten in Python by Yolian Amaro No newline at end of file
5 # Rewritten in python by Yolian Amaro No newline at end of file
6 # Python version 2.7 No newline at end of file
6 # Python version 2.7
7 # May 15, 2014 No newline at end of file
7 # May 15, 2014
8 # Jicamarca Radio Observatory No newline at end of file
8 # Jicamarca Radio Observatory
9 #---------------------------------------------------------- No newline at end of file
9 #----------------------------------------------------------
10
10
No newline at end of file
11 import time No newline at end of file
11 import math
No newline at end of file
12 import numpy as np No newline at end of file
13 import matplotlib.pyplot as plt
12 import matplotlib.pyplot as plt
No newline at end of file
13 from scipy.optimize import root
14 from scipy import linalg
No newline at end of file
No newline at end of file
14 No newline at end of file
15 import time No newline at end of file
16 from y_hysell96 import* No newline at end of file
15 from y_hysell96 import*
17 from deb4_basis import * No newline at end of file
16 from deb4_basis import *
18 from modelf import * No newline at end of file
17 from modelf import *
18 from irls_dn2 import * No newline at end of file
19 #from scipy.optimize import fsolve
19 #from scipy.optimize import fsolve
No newline at end of file
20
20 from scipy.optimize import root
No newline at end of file
No newline at end of file
21
21 import pywt
No newline at end of file
No newline at end of file
22 #-------------------------------------------------------------------------------------------------
22 from irls_dn2 import *
No newline at end of file
No newline at end of file
23 # Set parameters
23 No newline at end of file
No newline at end of file
24 #------------------------------------------------------------------------------------------------- No newline at end of file
24 No newline at end of file
25
25 ## Calculate Forward Model No newline at end of file
26 ## Calculate Forward Model
26 lambda1 = 6.0
27 lambda1 = 6.0
No newline at end of file
28 k = 2*np.pi/lambda1
27 k = 2*math.pi/lambda1
No newline at end of file
No newline at end of file
29
28
No newline at end of file
No newline at end of file
30 ## Magnetic Declination No newline at end of file
29 ## Calculate Magnetic Declination
No newline at end of file
30
No newline at end of file
31 # [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this
No newline at end of file
32
No newline at end of file
33 # or calculate it with the above function No newline at end of file
34 dec = -1.24 No newline at end of file
31 dec = -1.24
35
32
No newline at end of file
33 ## Loads Jicamarca antenna positions
36 # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt() No newline at end of file
No newline at end of file
34 antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False)
No newline at end of file
35
No newline at end of file
36 ## rx and ry -- for plotting purposes No newline at end of file
37 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
37 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
38 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
38 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
39
39
No newline at end of file
40 ## Plot of antenna positions No newline at end of file
40 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
41 [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
No newline at end of file
42 No newline at end of file
43 plt.figure(1) No newline at end of file
41 plt.figure(1)
44 plt.plot(rx, ry, 'ro') No newline at end of file
42 plt.plot(rx, ry, 'ro')
45 plt.draw() No newline at end of file
43 plt.draw()
46
44
No newline at end of file
45 ## Jicamarca is nominally at a 45 degree angle No newline at end of file
47 # Jicamarca is nominally at a 45 degree angle No newline at end of file
48 theta = 45 - dec; No newline at end of file
46 theta = 45 - dec;
49 No newline at end of file
47
48 ## Rotation matrix from antenna coord to magnetic coord (East North)
No newline at end of file
49 theta_rad = np.radians(theta) # trig functions take radians as argument
No newline at end of file
50 val1 = float( np.cos(theta_rad) )
No newline at end of file
51 val2 = float( np.sin(theta_rad) )
No newline at end of file
52 val3 = float( -1*np.sin(theta_rad))
No newline at end of file
53 val4 = float( np.cos(theta_rad) )
No newline at end of file
54 No newline at end of file
50 # Rotation matrix from antenna coord to magnetic coord (East North)
55 # Rotation matrix from antenna coord to magnetic coord (East North)
No newline at end of file
56 R = np.array( [[val1, val3], [val2, val4]] ) No newline at end of file
51 theta_rad = math.radians(theta) # trig functions take radians as argument
No newline at end of file
52 val1 = float( math.cos(theta_rad) )
No newline at end of file
53 val2 = float( math.sin(theta_rad) )
No newline at end of file
54 val3 = float( -1*math.sin(theta_rad))
No newline at end of file
55 val4 = float( math.cos(theta_rad) )
No newline at end of file
56
No newline at end of file
57 # Rotation matrix from antenna coord to magnetic coord (East North)
No newline at end of file
58 R = np.array( [[val1, val3], [val2, val4]] ); No newline at end of file
59 No newline at end of file
57
60 # Rotate antenna positions to magnetic coord.
58 # Rotate antenna positions to magnetic coord.
No newline at end of file
59 AR = np.dot(R.T, antpos) No newline at end of file
61 AR = np.dot(R.T, antpos); No newline at end of file
62 No newline at end of file
60
63 # Only take the East component No newline at end of file
61 # Only take the East component
64 r = AR[0,:]
62 r = AR[0,:]
No newline at end of file
63 r.sort() No newline at end of file
65 r.sort() # ROW VECTOR? No newline at end of file
66 No newline at end of file
64
67 # Truth model (high and low resolution)
65 # Truth model (high and low resolution)
No newline at end of file
66 Nt = (1024.0)*(16.0) # number of pixels in truth image: high resolution
68 Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
No newline at end of file
No newline at end of file
67 thbound = 9.0/180*np.pi # the width of the domain in angle space
69 thbound = 9.0/180*math.pi; # the width of the domain in angle space
No newline at end of file
No newline at end of file
68 thetat = np.linspace(-thbound, thbound,Nt) # image domain
70 thetat = np.linspace(-thbound, thbound,Nt) # image domain
No newline at end of file
No newline at end of file
69 thetat = thetat.T # transpose
71 thetat = np.transpose(thetat) # transpose # FUNCIONA??????????????????????????????
No newline at end of file
No newline at end of file
70 Nr = (256.0) # number of pixels in reconstructed image: low res No newline at end of file
72 Nr = (256.0); # number of pixels in reconstructed image: low res No newline at end of file
73 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
71 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
No newline at end of file
72 thetar = thetar.T # transpose
74 thetar = np.transpose(thetar) #transpose # FUNCIONA?????????????????????????????
No newline at end of file
No newline at end of file
73
75
No newline at end of file
No newline at end of file
74
76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
No newline at end of file
No newline at end of file
75 #-------------------------------------------------------------------------------------------------
77 # background constant b. No newline at end of file
No newline at end of file
76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b.
No newline at end of file
77 #------------------------------------------------------------------------------------------------- No newline at end of file
78 No newline at end of file
78
79 # Triple Gaussian No newline at end of file
79 # Triple Gaussian
80 # a = np.array([3, 5, 2]);
80 # a = np.array([3, 5, 2]);
No newline at end of file
81 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi, 7.0/180*np.pi]);
81 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
No newline at end of file
No newline at end of file
82 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi, 0.3/180*np.pi]);
82 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
No newline at end of file
No newline at end of file
83 # b = 0; # background No newline at end of file
83 # b = 0; # background No newline at end of file
84 No newline at end of file
84
85 # Double Gaussian No newline at end of file
85 # Double Gaussian
86 # a = np.array([3, 5]);
86 # a = np.array([3, 5]);
No newline at end of file
87 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi]);
87 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
No newline at end of file
No newline at end of file
88 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi]);
88 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
No newline at end of file
No newline at end of file
89 # b = 0; # background No newline at end of file
89 # b = 0; # background No newline at end of file
90 No newline at end of file
90
91 # Single Gaussian
91 # Single Gaussian
No newline at end of file
92 a = np.array( [3] )
92 a = np.array( [3] );
No newline at end of file
No newline at end of file
93 mu = np.array( [-3.0/180*np.pi] )
93 mu = np.array( [-3.0/180*math.pi] )
No newline at end of file
No newline at end of file
94 sig = np.array( [2.0/180*np.pi] )
94 sig = np.array( [2.0/180*math.pi] )
No newline at end of file
No newline at end of file
95 b = 0
95 b = 0;
No newline at end of file
No newline at end of file
96
96
No newline at end of file
No newline at end of file
97 # Empty matrices for factors
97 fact = np.zeros(shape=(Nt,1));
No newline at end of file
No newline at end of file
98 fact = np.zeros(shape=(Nt,1))
98 factr = np.zeros(shape=(Nr,1));
No newline at end of file
No newline at end of file
99 factr = np.zeros(shape=(Nr,1))
99 No newline at end of file
No newline at end of file
100
No newline at end of file
101 # DFT Kernels No newline at end of file
100 for i in range(0, a.size): No newline at end of file
102 for i in range(0, a.size):
101 temp = (-(thetat-mu[i])**2/(sig[i]**2)) No newline at end of file
103 temp = (-(thetat-mu[i])**2/(sig[i]**2))
102 tempr = (-(thetar-mu[i])**2/(sig[i]**2)) No newline at end of file
104 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
103 for j in range(0, temp.size):
105 for j in range(0, temp.size):
No newline at end of file
106 fact[j] = fact[j] + a[i]*np.exp(temp[j]); No newline at end of file
104 fact[j] = fact[j] + a[i]*math.exp(temp[j]); No newline at end of file
105 for m in range(0, tempr.size):
107 for m in range(0, tempr.size):
No newline at end of file
108 factr[m] = factr[m] + a[i]*np.exp(tempr[m]);
106 factr[m] = factr[m] + a[i]*math.exp(tempr[m]); No newline at end of file
No newline at end of file
109 No newline at end of file
107 fact = fact + b; No newline at end of file
110 fact = fact + b;
108 factr = factr + b; No newline at end of file
111 factr = factr + b;
109
112
No newline at end of file
113 # #-------------------------------------------------------------------------------------------------
110 # # model for f: Square pulse No newline at end of file
No newline at end of file
114 # # Model for f: Square pulse
No newline at end of file
115 # #------------------------------------------------------------------------------------------------- No newline at end of file
111 # for j in range(0, fact.size):
116 # for j in range(0, fact.size):
No newline at end of file
117 # if (theta > -5.0/180*np.pi and theta < 2.0/180*np.pi): No newline at end of file
112 # if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi): No newline at end of file
113 # fact[j] = 0 No newline at end of file
118 # fact[j] = 0
114 # else: No newline at end of file
119 # else:
115 # fact[j] = 1 No newline at end of file
120 # fact[j] = 1
116 # for k in range(0, factr.size):
121 # for k in range(0, factr.size):
No newline at end of file
122 # if (thetar[k] > -5.0/180*np.pi and thetar[k] < 2/180*np.pi):
117 # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
No newline at end of file
No newline at end of file
123 # factr[k] = 0 No newline at end of file
118 # fact[k] = 0 No newline at end of file
119 # else:
124 # else:
No newline at end of file
125 # factr[k] = 1
120 # fact[k] = 1
No newline at end of file
No newline at end of file
126
121 #
No newline at end of file
No newline at end of file
127 # #-------------------------------------------------------------------------------------------------
122 #
No newline at end of file
No newline at end of file
128 # # Model for f: Triangle pulse
123 # # model for f: triangle pulse
No newline at end of file
No newline at end of file
129 # #-------------------------------------------------------------------------------------------------
124 # mu = -1.0/180*math.pi;
No newline at end of file
No newline at end of file
130 # mu = -1.0/180*np.pi;
125 # sig = 5.0/180*math.pi; No newline at end of file
No newline at end of file
131 # sig = 5.0/180*np.pi; No newline at end of file
126 # wind1 = theta > mu-sig and theta < mu; No newline at end of file
132 # wind1 = theta > mu-sig and theta < mu;
127 # wind2 = theta < mu+sig and theta > mu; No newline at end of file
133 # wind2 = theta < mu+sig and theta > mu;
128 # fact = wind1 * (theta - (mu - sig)); No newline at end of file
134 # fact = wind1 * (theta - (mu - sig));
Show file before | Show file after | Hide comments
@@ -50,9 +50,11
50 lo = upfirdn(x, af[:,0], 1, 2); No newline at end of file
50 lo = upfirdn(x, af[:,0], 1, 2);
51 No newline at end of file
51
52
52
No newline at end of file
53 # # VERIFY THIS!!!!!!!!!!!!
53 # VERIFY THIS!!!!!!!!!!!!
No newline at end of file
No newline at end of file
54 # for i in range(0, L):
54 for i in range(0, L):
No newline at end of file
No newline at end of file
55 # lo[i] = lo[N/2+i] + lo[i];
55 lo[i] = lo[N/2+i] + lo[i]; No newline at end of file
No newline at end of file
56
No newline at end of file
57 lo[0:L-1] = lo[N/2+np.arange(0,L-1)] + lo[0:L-1] No newline at end of file
56 No newline at end of file
58
57 lo = lo[0:N/2]; No newline at end of file
59 lo = lo[0:N/2];
58 No newline at end of file
60
Show file before | Show file after | Hide comments
@@ -10,6 +10,7
10 from dualtree import * No newline at end of file
10 from dualtree import *
11 from idualtree import * No newline at end of file
11 from idualtree import *
12 No newline at end of file
12
13 # Debauchie 4 Wavelet No newline at end of file
13 def deb4_basis(N): No newline at end of file
14 def deb4_basis(N):
14 No newline at end of file
15
15 Psi = np.zeros(shape=(N,2*N+1)); No newline at end of file
16 Psi = np.zeros(shape=(N,2*N+1));
Show file before | Show file after | Hide comments
@@ -45,7 +45,7
45 #----------------------------------------# No newline at end of file
45 #----------------------------------------#
46 No newline at end of file
46
47 # normalization
47 # normalization
No newline at end of file
48 x = x/np.sqrt(2.0); No newline at end of file
48 x = x/np.sqrt(2); No newline at end of file
49 No newline at end of file
49
50 No newline at end of file
50
51 w = np.zeros(shape=(J+1), dtype=object) No newline at end of file
51 w = np.zeros(shape=(J+1), dtype=object)
Show file before | Show file after | Hide comments
@@ -41,6 +41,6
41 y2 = sfb(y2, w[0][1], Fsf[0,1]); No newline at end of file
41 y2 = sfb(y2, w[0][1], Fsf[0,1]);
42 No newline at end of file
42
43 # normalization
43 # normalization
No newline at end of file
44 y = (y1 + y2)/np.sqrt(2.0); No newline at end of file
44 y = (y1 + y2)/np.sqrt(2); No newline at end of file
45 No newline at end of file
45
46 return y No newline at end of file
46 return y
Show file before | Show file after | Hide comments
@@ -8,6 +8,7
8 from scipy import linalg No newline at end of file
8 from scipy import linalg
9 import scipy.sparse as sps No newline at end of file
9 import scipy.sparse as sps
10 import numpy as np No newline at end of file
10 import numpy as np
11 from numpy.linalg import norm No newline at end of file
11 No newline at end of file
12
12 def irls_dn(A,b,p,lambda1): No newline at end of file
13 def irls_dn(A,b,p,lambda1):
13 No newline at end of file
14
Show file before | Show file after | Hide comments
@@ -7,7 +7,11
7 from irls_dn import * No newline at end of file
7 from irls_dn import *
8 from scipy.optimize import fsolve No newline at end of file
8 from scipy.optimize import fsolve
9 import numpy as np
9 import numpy as np
No newline at end of file
10 from scipy.optimize import *
10 from scipy.optimize import root No newline at end of file
No newline at end of file
11 from dogleg import *
No newline at end of file
12 from numpy.linalg import norm
No newline at end of file
13
No newline at end of file
14 import matplotlib.pyplot as plt No newline at end of file
11 No newline at end of file
15
12 def irls_dn2(A,b,p,G): No newline at end of file
16 def irls_dn2(A,b,p,G):
13 No newline at end of file
17
Show file before | Show file after | Hide comments
@@ -45,8 +45,12
45 No newline at end of file
45
46 lo = upfirdn(lo, sf[:,0], 2, 1); No newline at end of file
46 lo = upfirdn(lo, sf[:,0], 2, 1);
47 hi = upfirdn(hi, sf[:,1], 2, 1); No newline at end of file
47 hi = upfirdn(hi, sf[:,1], 2, 1);
48
No newline at end of file
49 lo = lo[0:lo.size-1]
No newline at end of file
50 hi = hi[0:hi.size-1]
No newline at end of file
51 No newline at end of file
48 y = lo + hi;
52 y = lo + hi;
No newline at end of file
53 y[0:L-2] = y[0:L-2] + y[N+ np.arange(0,L-2)]; #CHECK IF ARANGE IS CORRECT No newline at end of file
49 y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT No newline at end of file
50 y = y[0:N]; No newline at end of file
54 y = y[0:N];
51 No newline at end of file
55
52 #print 'y en sbf\n', y.shape No newline at end of file
56 #print 'y en sbf\n', y.shape
General Comments 0
You need to be logged in to leave comments. Login now