| @@ -0,0 +1,46 | |||
|
|
1 | ''' | |
|
|
No newline at end of file | ||
|
|
2 | Created on Jun 5, 2014 | |
|
|
No newline at end of file | ||
|
|
3 | ||
|
|
No newline at end of file | ||
|
|
4 | @author: Yolian Amaro | |
|
|
No newline at end of file | ||
|
|
5 | ''' | |
|
|
No newline at end of file | ||
|
|
6 | ||
|
|
No newline at end of file | ||
|
|
7 | from sfb import * | |
|
|
No newline at end of file | ||
|
|
8 | ||
|
|
No newline at end of file | ||
|
|
9 | def idualtree(w, J, Fsf, sf): | |
|
|
No newline at end of file | ||
|
|
10 | ||
|
|
No newline at end of file | ||
|
|
11 | # Inverse Dual-tree Complex DWT | |
|
|
No newline at end of file | ||
|
|
12 | # | |
|
|
No newline at end of file | ||
|
|
13 | # USAGE: | |
|
|
No newline at end of file | ||
|
|
14 | # y = idualtree(w, J, Fsf, sf) | |
|
|
No newline at end of file | ||
|
|
15 | # INPUT: | |
|
|
No newline at end of file | ||
|
|
16 | # w - DWT coefficients | |
|
|
No newline at end of file | ||
|
|
17 | # J - number of stages | |
|
|
No newline at end of file | ||
|
|
18 | # Fsf - synthesis filters for the last stage | |
|
|
No newline at end of file | ||
|
|
19 | # sf - synthesis filters for preceeding stages | |
|
|
No newline at end of file | ||
|
|
20 | # OUTUT: | |
|
|
No newline at end of file | ||
|
|
21 | # y - output signal | |
|
|
No newline at end of file | ||
|
|
22 | # See dualtree | |
|
|
No newline at end of file | ||
|
|
23 | # | |
|
|
No newline at end of file | ||
|
|
24 | # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY | |
|
|
No newline at end of file | ||
|
|
25 | # http://taco.poly.edu/WaveletSoftware/ | |
|
|
No newline at end of file | ||
|
|
26 | ||
|
|
No newline at end of file | ||
|
|
27 | # Tree 1 | |
|
|
No newline at end of file | ||
|
|
28 | y1 = w[J][0]; | |
|
|
No newline at end of file | ||
|
|
29 | ||
|
|
No newline at end of file | ||
|
|
30 | for j in range (J-1, 0, -1): | |
|
|
No newline at end of file | ||
|
|
31 | y1 = sfb(y1, w[j][0], sf[0,0]); | |
|
|
No newline at end of file | ||
|
|
32 | ||
|
|
No newline at end of file | ||
|
|
33 | y1 = sfb(y1, w[0][0], Fsf[0,0]); | |
|
|
No newline at end of file | ||
|
|
34 | ||
|
|
No newline at end of file | ||
|
|
35 | # Tree 2 | |
|
|
No newline at end of file | ||
|
|
36 | y2 = w[J][1]; | |
|
|
No newline at end of file | ||
|
|
37 | ||
|
|
No newline at end of file | ||
|
|
38 | for j in range (J-1, 0, -1): | |
|
|
No newline at end of file | ||
|
|
39 | y2 = sfb(y2, w[j][2], sf[0,1]); | |
|
|
No newline at end of file | ||
|
|
40 | ||
|
|
No newline at end of file | ||
|
|
41 | y2 = sfb(y2, w[0][1], Fsf[0,1]); | |
|
|
No newline at end of file | ||
|
|
42 | ||
|
|
No newline at end of file | ||
|
|
43 | # normalization | |
|
|
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 | ||
|
|
46 | return y No newline at end of file | |
| @@ -0,0 +1,68 | |||
|
|
1 | ''' | |
|
|
No newline at end of file | ||
|
|
2 | Created on Jun 5, 2014 | |
|
|
No newline at end of file | ||
|
|
3 | ||
|
|
No newline at end of file | ||
|
|
4 | @author: Yolian Amaro | |
|
|
No newline at end of file | ||
|
|
5 | ''' | |
|
|
No newline at end of file | ||
|
|
6 | ||
|
|
No newline at end of file | ||
|
|
7 | from multirate import * | |
|
|
No newline at end of file | ||
|
|
8 | import numpy as np | |
|
|
No newline at end of file | ||
|
|
9 | from cshift import * | |
|
|
No newline at end of file | ||
|
|
10 | ||
|
|
No newline at end of file | ||
|
|
11 | def sfb(lo, hi, sf): | |
|
|
No newline at end of file | ||
|
|
12 | ||
|
|
No newline at end of file | ||
|
|
13 | # Synthesis filter bank | |
|
|
No newline at end of file | ||
|
|
14 | # | |
|
|
No newline at end of file | ||
|
|
15 | # USAGE: | |
|
|
No newline at end of file | ||
|
|
16 | # y = sfb(lo, hi, sf) | |
|
|
No newline at end of file | ||
|
|
17 | # INPUT: | |
|
|
No newline at end of file | ||
|
|
18 | # lo - low frqeuency input | |
|
|
No newline at end of file | ||
|
|
19 | # hi - high frequency input | |
|
|
No newline at end of file | ||
|
|
20 | # sf - synthesis filters | |
|
|
No newline at end of file | ||
|
|
21 | # sf(:, 1) - lowpass filter (even length) | |
|
|
No newline at end of file | ||
|
|
22 | # sf(:, 2) - highpass filter (even length) | |
|
|
No newline at end of file | ||
|
|
23 | # OUTPUT: | |
|
|
No newline at end of file | ||
|
|
24 | # y - output signal | |
|
|
No newline at end of file | ||
|
|
25 | # See also afb | |
|
|
No newline at end of file | ||
|
|
26 | # | |
|
|
No newline at end of file | ||
|
|
27 | # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY | |
|
|
No newline at end of file | ||
|
|
28 | # http://taco.poly.edu/WaveletSoftware/ | |
|
|
No newline at end of file | ||
|
|
29 | ||
|
|
No newline at end of file | ||
|
|
30 | N = 2*lo.size; | |
|
|
No newline at end of file | ||
|
|
31 | L = sf.size/2; | |
|
|
No newline at end of file | ||
|
|
32 | #print 'N', N | |
|
|
No newline at end of file | ||
|
|
33 | #print 'sf', sf | |
|
|
No newline at end of file | ||
|
|
34 | ||
|
|
No newline at end of file | ||
|
|
35 | ||
|
|
No newline at end of file | ||
|
|
36 | #print 'sf[:,0]', sf[:,0].shape | |
|
|
No newline at end of file | ||
|
|
37 | #print 'sf[:,1]', sf[:,1].shape | |
|
|
No newline at end of file | ||
|
|
38 | #print 'sbf hi', hi.shape | |
|
|
No newline at end of file | ||
|
|
39 | ||
|
|
No newline at end of file | ||
|
|
40 | ||
|
|
No newline at end of file | ||
|
|
41 | ||
|
|
No newline at end of file | ||
|
|
42 | # Need to change format for upfirdn funct: | |
|
|
No newline at end of file | ||
|
|
43 | lo = lo.T.conj() | |
|
|
No newline at end of file | ||
|
|
44 | lo = lo.reshape(lo.size) | |
|
|
No newline at end of file | ||
|
|
45 | ||
|
|
No newline at end of file | ||
|
|
46 | print 'sfb hi', hi | |
|
|
No newline at end of file | ||
|
|
47 | ||
|
|
No newline at end of file | ||
|
|
48 | # Need to change format for upfirdn funct: | |
|
|
No newline at end of file | ||
|
|
49 | hi = hi.T.conj() | |
|
|
No newline at end of file | ||
|
|
50 | hi = hi.reshape(hi.size) | |
|
|
No newline at end of file | ||
|
|
51 | ||
|
|
No newline at end of file | ||
|
|
52 | #hi = hi.reshape(1, hi.size) | |
|
|
No newline at end of file | ||
|
|
53 | ||
|
|
No newline at end of file | ||
|
|
54 | lo = upfirdn(lo, sf[:,0], 2, 1); | |
|
|
No newline at end of file | ||
|
|
55 | hi = upfirdn(hi, sf[:,1], 2, 1); | |
|
|
No newline at end of file | ||
|
|
56 | y = lo + hi; | |
|
|
No newline at end of file | ||
|
|
57 | 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 | ||
|
|
58 | y = y[0:N]; | |
|
|
No newline at end of file | ||
|
|
59 | ||
|
|
No newline at end of file | ||
|
|
60 | print 'y en sbf\n', y.shape | |
|
|
No newline at end of file | ||
|
|
61 | ||
|
|
No newline at end of file | ||
|
|
62 | y = y.reshape(1, y.size) | |
|
|
No newline at end of file | ||
|
|
63 | print 'y en sbf\n', y.shape | |
|
|
No newline at end of file | ||
|
|
64 | ||
|
|
No newline at end of file | ||
|
|
65 | y = cshift(y, 1-L/2); | |
|
|
No newline at end of file | ||
|
|
66 | ||
|
|
No newline at end of file | ||
|
|
67 | return y; | |
|
|
No newline at end of file | ||
|
|
68 | No newline at end of file | |
| @@ -1,65 +1,63 | |||
|
|
1 | 1 | ''' No newline at end of file |
|
|
2 | 2 | Created on May 26, 2014 No newline at end of file |
|
|
3 | 3 | No newline at end of file |
|
|
4 | 4 | @author: Yolian Amaro No newline at end of file |
|
|
5 | 5 | ''' No newline at end of file |
|
|
6 | 6 | No newline at end of file |
|
|
7 | 7 | import pywt No newline at end of file |
|
|
8 | 8 | import numpy as np No newline at end of file |
|
|
9 | 9 | No newline at end of file |
|
|
10 | 10 | def FSfarras(): No newline at end of file |
|
|
11 | 11 | #function [af, sf] = FSfarras No newline at end of file |
|
|
12 | 12 | No newline at end of file |
|
|
13 | 13 | # Farras filters organized for the dual-tree No newline at end of file |
|
|
14 | 14 | # complex DWT. No newline at end of file |
|
|
15 | 15 | # No newline at end of file |
|
|
16 | 16 | # USAGE: No newline at end of file |
|
|
17 | 17 | # [af, sf] = FSfarras No newline at end of file |
|
|
18 | 18 | # OUTPUT: No newline at end of file |
|
|
19 | 19 | # af{i}, i = 1,2 - analysis filters for tree i No newline at end of file |
|
|
20 | 20 | # sf{i}, i = 1,2 - synthesis filters for tree i No newline at end of file |
|
|
21 | 21 | # See farras, dualtree, dualfilt1. No newline at end of file |
|
|
22 | 22 | # No newline at end of file |
|
|
23 | 23 | # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY No newline at end of file |
|
|
24 | 24 | # http://taco.poly.edu/WaveletSoftware/ No newline at end of file |
|
|
25 | 25 | # No newline at end of file |
|
|
26 | 26 | # Translated to Python by Yolian Amaro No newline at end of file |
|
|
27 | 27 | No newline at end of file |
|
|
28 | 28 | No newline at end of file |
|
|
29 | 29 | No newline at end of file |
|
|
30 | 30 | a1 = np.array( [ No newline at end of file |
|
|
31 | 31 | [ 0, 0], No newline at end of file |
|
|
32 | 32 |
[ |
|
|
33 | 33 | [ 0.08838834764832, 0.01122679215254], No newline at end of file |
|
|
34 | 34 |
[ 0.69587998903400, |
|
|
35 | 35 | [ 0.69587998903400, 0.08838834764832], No newline at end of file |
|
|
36 | 36 |
[ |
|
|
37 | 37 | [-0.08838834764832, 0.69587998903400], No newline at end of file |
|
|
38 | 38 | [ 0.01122679215254, -0.08838834764832], No newline at end of file |
|
|
39 | 39 | [ 0.01122679215254, -0.08838834764832], No newline at end of file |
|
|
40 | 40 | [0, 0] No newline at end of file |
|
|
41 | 41 |
|
|
|
42 | 42 | No newline at end of file |
|
|
43 | 43 | a2 = np.array([ No newline at end of file |
|
|
44 | 44 | [ 0.01122679215254, 0], No newline at end of file |
|
|
45 | 45 | [ 0.01122679215254, 0], No newline at end of file |
|
|
46 | 46 |
[ |
|
|
47 | 47 | [ 0.08838834764832, -0.08838834764832], No newline at end of file |
|
|
48 | 48 |
[ 0.69587998903400, |
|
|
49 | 49 | [ 0.69587998903400, -0.69587998903400], No newline at end of file |
|
|
50 | 50 |
[ |
|
|
51 | 51 | [-0.08838834764832, 0.08838834764832], No newline at end of file |
|
|
52 | 52 |
[ 0, |
|
|
53 | 53 | [ 0, -0.01122679215254] No newline at end of file |
|
|
54 | 54 |
|
|
|
55 | 55 |
|
|
|
No newline at end of file | ||
|
|
56 | #print a2.shape No newline at end of file | |
|
|
57 | 56 | No newline at end of file |
|
|
58 | 57 | af = np.array([ [a1,a2] ], dtype=object) No newline at end of file |
|
|
59 | 58 | No newline at end of file |
|
|
60 | 59 |
s |
|
|
61 | 60 | s2 = a2[::-1] No newline at end of file |
|
|
62 | 61 | No newline at end of file |
|
|
63 | 62 | sf = np.array([ [s1,s2] ], dtype=object) No newline at end of file |
|
|
64 | 63 | No newline at end of file |
|
|
65 | 64 | return af, sf No newline at end of file |
| @@ -1,468 +1,471 | |||
|
|
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 numpy as np No newline at end of file |
|
|
13 | 13 | import matplotlib.pyplot as plt No newline at end of file |
|
|
14 | 14 | from scipy import linalg No newline at end of file |
|
|
15 | 15 | import time No newline at end of file |
|
|
16 | 16 | from y_hysell96 import* No newline at end of file |
|
|
17 | 17 | from deb4_basis import * No newline at end of file |
|
|
18 | 18 | from modelf import * No newline at end of file |
|
|
19 | 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 |
|
|
21 | 21 | import pywt No newline at end of file |
|
|
22 | 22 | from irls_dn2 import * No newline at end of file |
|
|
23 | 23 | No newline at end of file |
|
|
24 | 24 | No newline at end of file |
|
|
25 | 25 | ## Calculate Forward Model No newline at end of file |
|
|
26 | 26 | lambda1 = 6.0 No newline at end of file |
|
|
27 | 27 | k = 2*math.pi/lambda1 No newline at end of file |
|
|
28 | 28 | No newline at end of file |
|
|
29 | 29 | ## Calculate Magnetic Declination No newline at end of file |
|
|
30 | 30 | No newline at end of file |
|
|
31 | 31 | # [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this No newline at end of file |
|
|
32 | 32 | No newline at end of file |
|
|
33 | 33 | # or calculate it with the above function No newline at end of file |
|
|
34 | 34 | dec = -1.24 No newline at end of file |
|
|
35 | 35 | No newline at end of file |
|
|
36 | 36 | # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt() No newline at end of file |
|
|
37 | 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 |
|
|
38 | 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 |
|
|
39 | 39 | No newline at end of file |
|
|
40 | 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 | 41 | [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] ) No newline at end of file |
|
|
42 | 42 | No newline at end of file |
|
|
43 | 43 | plt.figure(1) No newline at end of file |
|
|
44 | 44 | plt.plot(rx, ry, 'ro') No newline at end of file |
|
|
45 | 45 | plt.draw() No newline at end of file |
|
|
46 | 46 | No newline at end of file |
|
|
47 | 47 | # Jicamarca is nominally at a 45 degree angle No newline at end of file |
|
|
48 | 48 | theta = 45 - dec; No newline at end of file |
|
|
49 | 49 | No newline at end of file |
|
|
50 | 50 | # Rotation matrix from antenna coord to magnetic coord (East North) No newline at end of file |
|
|
51 | 51 | theta_rad = math.radians(theta) # trig functions take radians as argument No newline at end of file |
|
|
52 | 52 | val1 = float( math.cos(theta_rad) ) No newline at end of file |
|
|
53 | 53 | val2 = float( math.sin(theta_rad) ) No newline at end of file |
|
|
54 | 54 | val3 = float( -1*math.sin(theta_rad)) No newline at end of file |
|
|
55 | 55 | val4 = float( math.cos(theta_rad) ) No newline at end of file |
|
|
56 | 56 | No newline at end of file |
|
|
57 | 57 | # Rotation matrix from antenna coord to magnetic coord (East North) No newline at end of file |
|
|
58 | 58 | R = np.array( [[val1, val3], [val2, val4]] ); No newline at end of file |
|
|
59 | 59 | No newline at end of file |
|
|
60 | 60 | # Rotate antenna positions to magnetic coord. No newline at end of file |
|
|
61 | 61 | AR = np.dot(R.T, antpos); No newline at end of file |
|
|
62 | 62 | No newline at end of file |
|
|
63 | 63 | # Only take the East component No newline at end of file |
|
|
64 | 64 | r = AR[0,:] No newline at end of file |
|
|
65 | 65 | r.sort() # ROW VECTOR? No newline at end of file |
|
|
66 | 66 | No newline at end of file |
|
|
67 | 67 | # Truth model (high and low resolution) No newline at end of file |
|
|
68 | 68 | Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution No newline at end of file |
|
|
69 | 69 | thbound = 9.0/180*math.pi; # the width of the domain in angle space No newline at end of file |
|
|
70 | 70 | thetat = np.linspace(-thbound, thbound,Nt) # image domain No newline at end of file |
|
|
71 | 71 | thetat = np.transpose(thetat) # transpose # FUNCIONA?????????????????????????????? No newline at end of file |
|
|
72 | 72 | Nr = (256.0); # number of pixels in reconstructed image: low res No newline at end of file |
|
|
73 | 73 | thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain No newline at end of file |
|
|
74 | 74 | thetar = np.transpose(thetar) #transpose # FUNCIONA????????????????????????????? No newline at end of file |
|
|
75 | 75 | No newline at end of file |
|
|
76 | 76 | # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and No newline at end of file |
|
|
77 | 77 | # background constant b. No newline at end of file |
|
|
78 | 78 | No newline at end of file |
|
|
79 | 79 | # Triple Gaussian No newline at end of file |
|
|
80 | 80 | # a = np.array([3, 5, 2]); No newline at end of file |
|
|
81 | 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 |
|
|
82 | 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 |
|
|
83 | 83 | # b = 0; # background No newline at end of file |
|
|
84 | 84 | No newline at end of file |
|
|
85 | 85 | # Double Gaussian No newline at end of file |
|
|
86 | 86 | # a = np.array([3, 5]); No newline at end of file |
|
|
87 | 87 | # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]); No newline at end of file |
|
|
88 | 88 | # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]); No newline at end of file |
|
|
89 | 89 | # b = 0; # background No newline at end of file |
|
|
90 | 90 | No newline at end of file |
|
|
91 | 91 | # Single Gaussian No newline at end of file |
|
|
92 | 92 | a = np.array( [3] ); No newline at end of file |
|
|
93 | 93 | mu = np.array( [-3.0/180*math.pi] ) No newline at end of file |
|
|
94 | 94 | sig = np.array( [2.0/180*math.pi] ) No newline at end of file |
|
|
95 | 95 | b = 0; No newline at end of file |
|
|
96 | 96 | No newline at end of file |
|
|
97 | 97 | fact = np.zeros(shape=(Nt,1)); No newline at end of file |
|
|
98 | 98 | factr = np.zeros(shape=(Nr,1)); No newline at end of file |
|
|
99 | 99 | No newline at end of file |
|
|
100 | 100 | for i in range(0, a.size): No newline at end of file |
|
|
101 | 101 | temp = (-(thetat-mu[i])**2/(sig[i]**2)) No newline at end of file |
|
|
102 | 102 | tempr = (-(thetar-mu[i])**2/(sig[i]**2)) No newline at end of file |
|
|
103 | 103 | for j in range(0, temp.size): No newline at end of file |
|
|
104 | 104 | fact[j] = fact[j] + a[i]*math.exp(temp[j]); No newline at end of file |
|
|
105 | 105 | for m in range(0, tempr.size): No newline at end of file |
|
|
106 | 106 | factr[m] = factr[m] + a[i]*math.exp(tempr[m]); No newline at end of file |
|
|
107 | 107 | fact = fact + b; No newline at end of file |
|
|
108 | 108 | factr = factr + b; No newline at end of file |
|
|
109 | 109 | No newline at end of file |
|
|
110 | 110 | # # model for f: Square pulse No newline at end of file |
|
|
111 | 111 | # for j in range(0, fact.size): No newline at end of file |
|
|
112 | 112 | # if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi): No newline at end of file |
|
|
113 | 113 | # fact[j] = 0 No newline at end of file |
|
|
114 | 114 | # else: No newline at end of file |
|
|
115 | 115 | # fact[j] = 1 No newline at end of file |
|
|
116 | 116 | # for k in range(0, factr.size): No newline at end of file |
|
|
117 | 117 | # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi): No newline at end of file |
|
|
118 | 118 | # fact[k] = 0 No newline at end of file |
|
|
119 | 119 | # else: No newline at end of file |
|
|
120 | 120 | # fact[k] = 1 No newline at end of file |
|
|
121 | 121 | # No newline at end of file |
|
|
122 | 122 | # No newline at end of file |
|
|
123 | 123 | # # model for f: triangle pulse No newline at end of file |
|
|
124 | 124 | # mu = -1.0/180*math.pi; No newline at end of file |
|
|
125 | 125 | # sig = 5.0/180*math.pi; No newline at end of file |
|
|
126 | 126 | # wind1 = theta > mu-sig and theta < mu; No newline at end of file |
|
|
127 | 127 | # wind2 = theta < mu+sig and theta > mu; No newline at end of file |
|
|
128 | 128 | # fact = wind1 * (theta - (mu - sig)); No newline at end of file |
|
|
129 | 129 | # factr = wind1 * (thetar - (mu - sig)); No newline at end of file |
|
|
130 | 130 | # fact = fact + wind2 * (-(theta-(mu+sig))); No newline at end of file |
|
|
131 | 131 | # factr = factr + wind2 * (-(thetar-(mu+sig))); No newline at end of file |
|
|
132 | 132 | No newline at end of file |
|
|
133 | 133 | No newline at end of file |
|
|
134 | 134 | # fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1 No newline at end of file |
|
|
135 | 135 | I = sum(fact)[0]; No newline at end of file |
|
|
136 | 136 | fact = fact/I; # normalize to sum(f)==1 No newline at end of file |
|
|
137 | 137 | factr = factr/I; # normalize to sum(f)==1 No newline at end of file |
|
|
138 | 138 | #plt.figure() No newline at end of file |
|
|
139 | 139 | #plt.plot(thetat,fact,'r'); No newline at end of file |
|
|
140 | 140 | #plt.plot(thetar,factr,'k.'); No newline at end of file |
|
|
141 | 141 | #xlim([min(thetat) max(thetat)]); No newline at end of file |
|
|
142 | 142 | No newline at end of file |
|
|
143 | 143 | #x = np.linspace(thetat.min(), thetat.max) ???? No newline at end of file |
|
|
144 | 144 | #for i in range(0, thetat.size): No newline at end of file |
|
|
145 | 145 | plt.figure(2) No newline at end of file |
|
|
146 | 146 | plt.plot(thetat, fact, 'r--') No newline at end of file |
|
|
147 | 147 | plt.plot(thetar, factr, 'ro') No newline at end of file |
|
|
148 | 148 | plt.draw() No newline at end of file |
|
|
149 | 149 | # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO No newline at end of file |
|
|
150 | 150 | No newline at end of file |
|
|
151 | 151 | No newline at end of file |
|
|
152 | 152 | ## No newline at end of file |
|
|
153 | 153 | # Control the type and number of inversions with: No newline at end of file |
|
|
154 | 154 | # SNRdBvec: the SNRs that will be used. No newline at end of file |
|
|
155 | 155 | # NN: the number of trials for each SNR No newline at end of file |
|
|
156 | 156 | No newline at end of file |
|
|
157 | 157 | #SNRdBvec = np.linspace(5,20,10); No newline at end of file |
|
|
158 | 158 | SNRdBvec = np.array([15]); No newline at end of file |
|
|
159 | 159 | NN = 1; # number of trial at each SNR No newline at end of file |
|
|
160 | 160 | No newline at end of file |
|
|
161 | 161 | # if using vector arguments should be: (4,SNRdBvec.size,NN) No newline at end of file |
|
|
162 | 162 | corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file |
|
|
163 | 163 | corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file |
|
|
164 | 164 | rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) No newline at end of file |
|
|
165 | 165 | No newline at end of file |
|
|
166 | 166 | for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector No newline at end of file |
|
|
167 | 167 | for Ni in range(0, NN): No newline at end of file |
|
|
168 | 168 | SNRdB = SNRdBvec[snri]; No newline at end of file |
|
|
169 | 169 | SNR = 10**(SNRdB/10.0); No newline at end of file |
|
|
170 | 170 | No newline at end of file |
|
|
171 | 171 | # Calculate cross-correlation matrix (Fourier components of image) No newline at end of file |
|
|
172 | 172 | # This is an inefficient way to do this. No newline at end of file |
|
|
173 | 173 | R = np.zeros(shape=(r.size, r.size), dtype=object); No newline at end of file |
|
|
174 | 174 | No newline at end of file |
|
|
175 | 175 | for i1 in range(0, r.size): No newline at end of file |
|
|
176 | 176 | for i2 in range(0,r.size): No newline at end of file |
|
|
177 | 177 | 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 |
|
|
178 | 178 | R[i1,i2] = sum(R[i1,i2]) No newline at end of file |
|
|
179 | 179 | No newline at end of file |
|
|
180 | 180 | # Add uncertainty No newline at end of file |
|
|
181 | 181 | # This is an ad-hoc way of adding "noise". It models some combination of No newline at end of file |
|
|
182 | 182 | # receiver noise and finite integration times. We could use a more No newline at end of file |
|
|
183 | 183 | # advanced model (like in Yu et al 2000) in the future. No newline at end of file |
|
|
184 | 184 | No newline at end of file |
|
|
185 | 185 | # This is a way of adding noise while maintaining the No newline at end of file |
|
|
186 | 186 | # positive-semi-definiteness of the matrix. No newline at end of file |
|
|
187 | 187 | No newline at end of file |
|
|
188 | 188 | U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R No newline at end of file |
|
|
189 | 189 | No newline at end of file |
|
|
190 | 190 | sigma_noise = (np.linalg.norm(U,'fro')/SNR); No newline at end of file |
|
|
191 | 191 | No newline at end of file |
|
|
192 | 192 | temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5) No newline at end of file |
|
|
193 | 193 | temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5) No newline at end of file |
|
|
194 | 194 | temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's No newline at end of file |
|
|
195 | 195 | temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0) No newline at end of file |
|
|
196 | 196 | No newline at end of file |
|
|
197 | 197 | nz = np.multiply(temp4, temp3) No newline at end of file |
|
|
198 | 198 | No newline at end of file |
|
|
199 | 199 | #---------------------- Eliminar esto:------------------------------------------ No newline at end of file |
|
|
200 | 200 | #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int)) No newline at end of file |
|
|
201 | 201 | #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int)) No newline at end of file |
|
|
202 | 202 | #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 )); No newline at end of file |
|
|
203 | 203 | #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 |
|
|
204 | 204 | #-------------------------------------------------------------------------------- No newline at end of file |
|
|
205 | 205 | No newline at end of file |
|
|
206 | 206 | Unz = U + nz; No newline at end of file |
|
|
207 | 207 | Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R No newline at end of file |
|
|
208 | 208 | plt.figure(3); No newline at end of file |
|
|
209 | 209 | plt.pcolor(abs(Rnz)); No newline at end of file |
|
|
210 | 210 | plt.colorbar(); No newline at end of file |
|
|
211 | 211 | No newline at end of file |
|
|
212 | 212 | # Fourier Inversion ################### No newline at end of file |
|
|
213 | 213 | f_fourier = np.zeros(shape=(Nr,1), dtype=complex); No newline at end of file |
|
|
214 | 214 | No newline at end of file |
|
|
215 | 215 | for i in range(0, thetar.size): No newline at end of file |
|
|
216 | 216 | th = thetar[i]; No newline at end of file |
|
|
217 | 217 | w = np.exp(1j*k*np.dot(r,np.sin(th))); No newline at end of file |
|
|
218 | 218 | No newline at end of file |
|
|
219 | 219 | temp = np.dot(w.T.conj(),U) No newline at end of file |
|
|
220 | 220 | No newline at end of file |
|
|
221 | 221 | f_fourier[i] = np.dot(temp, w); No newline at end of file |
|
|
222 | 222 | No newline at end of file |
|
|
223 | 223 | f_fourier = f_fourier.real; # get rid of numerical imaginary noise No newline at end of file |
|
|
224 | 224 | No newline at end of file |
|
|
225 | 225 | #print f_fourier No newline at end of file |
|
|
226 | 226 | No newline at end of file |
|
|
227 | 227 | No newline at end of file |
|
|
228 | 228 | # Capon Inversion ###################### No newline at end of file |
|
|
229 | 229 | No newline at end of file |
|
|
230 | 230 | f_capon = np.zeros(shape=(Nr,1)); No newline at end of file |
|
|
231 | 231 | No newline at end of file |
|
|
232 | 232 | tic_capon = time.time(); No newline at end of file |
|
|
233 | 233 | No newline at end of file |
|
|
234 | 234 | for i in range(0, thetar.size): No newline at end of file |
|
|
235 | 235 | th = thetar[i]; No newline at end of file |
|
|
236 | 236 | w = np.exp(1j*k*np.dot(r,np.sin(th))); No newline at end of file |
|
|
237 | 237 | f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real) No newline at end of file |
|
|
238 | 238 | No newline at end of file |
|
|
239 | 239 | No newline at end of file |
|
|
240 | 240 | toc_capon = time.time() No newline at end of file |
|
|
241 | 241 | No newline at end of file |
|
|
242 | 242 | elapsed_time_capon = toc_capon - tic_capon; No newline at end of file |
|
|
243 | 243 | No newline at end of file |
|
|
244 | 244 | f_capon = f_capon.real; # get rid of numerical imaginary noise No newline at end of file |
|
|
245 | 245 | No newline at end of file |
|
|
246 | 246 | # MaxEnt Inversion ##################### No newline at end of file |
|
|
247 | 247 | No newline at end of file |
|
|
248 | 248 | # create the appropriate sensing matrix (split into real and imaginary # parts) No newline at end of file |
|
|
249 | 249 | M = (r.size-1)*(r.size); No newline at end of file |
|
|
250 | 250 | Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix No newline at end of file |
|
|
251 | 251 | Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction No newline at end of file |
|
|
252 | 252 | No newline at end of file |
|
|
253 | 253 | # need to re-index our measurements from matrix R into vector g No newline at end of file |
|
|
254 | 254 | g = np.zeros(shape=(M,1)); No newline at end of file |
|
|
255 | 255 | gnz = np.zeros(shape=(M,1)); # noisy version of g No newline at end of file |
|
|
256 | 256 | No newline at end of file |
|
|
257 | 257 | # triangular indexing to perform this re-indexing No newline at end of file |
|
|
258 | 258 | T = np.ones(shape=(r.size,r.size)); No newline at end of file |
|
|
259 | 259 | [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing No newline at end of file |
|
|
260 | 260 | No newline at end of file |
|
|
261 | 261 | # build H No newline at end of file |
|
|
262 | 262 | for i1 in range(0, r.size): No newline at end of file |
|
|
263 | 263 | for i2 in range(i1+1, r.size): No newline at end of file |
|
|
264 | 264 | idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward No newline at end of file |
|
|
265 | 265 | idx1 = 2*idx; # because index starts at 0 No newline at end of file |
|
|
266 | 266 | idx2 = 2*idx+1; No newline at end of file |
|
|
267 | 267 | Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T; No newline at end of file |
|
|
268 | 268 | Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T; No newline at end of file |
|
|
269 | 269 | Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt; No newline at end of file |
|
|
270 | 270 | Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt; No newline at end of file |
|
|
271 | 271 | g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later No newline at end of file |
|
|
272 | 272 | g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again No newline at end of file |
|
|
273 | 273 | gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt; No newline at end of file |
|
|
274 | 274 | gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt; No newline at end of file |
|
|
275 | 275 | No newline at end of file |
|
|
276 | 276 | # inversion No newline at end of file |
|
|
277 | 277 | F = Nr/Nt; # normalization No newline at end of file |
|
|
278 | 278 | sigma = 1; # set to 1 because the difference is accounted for in G No newline at end of file |
|
|
279 | 279 | No newline at end of file |
|
|
280 | 280 | ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below No newline at end of file |
|
|
281 | 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 |
|
|
282 | 282 | No newline at end of file |
|
|
283 | 283 | tic_maxent = time.time(); No newline at end of file |
|
|
284 | 284 | No newline at end of file |
|
|
285 | 285 | lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything) No newline at end of file |
|
|
286 | 286 | No newline at end of file |
|
|
287 | 287 | toc_maxent = time.time() No newline at end of file |
|
|
288 | 288 | elapsed_time_maxent = toc_maxent - tic_maxent; No newline at end of file |
|
|
289 | 289 | No newline at end of file |
|
|
290 | 290 | # Whitened solution No newline at end of file |
|
|
291 | 291 | def myfun(lambda1): No newline at end of file |
|
|
292 | 292 | return y_hysell96(lambda1,gnz,sigma,F,G,Hr); No newline at end of file |
|
|
293 | 293 | No newline at end of file |
|
|
294 | 294 | tic_maxEnt = time.time(); No newline at end of file |
|
|
295 | 295 | No newline at end of file |
|
|
296 | 296 | #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000); No newline at end of file |
|
|
297 | 297 | lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14); No newline at end of file |
|
|
298 | 298 | No newline at end of file |
|
|
299 | 299 | #print lambda1 No newline at end of file |
|
|
300 | 300 | #print lambda1.x No newline at end of file |
|
|
301 | 301 | No newline at end of file |
|
|
302 | 302 | lambda1 = lambda1.x; No newline at end of file |
|
|
303 | 303 | No newline at end of file |
|
|
304 | 304 | toc_maxEnt = time.time(); No newline at end of file |
|
|
305 | 305 | f_maxent = modelf(lambda1, Hr, F); No newline at end of file |
|
|
306 | 306 | ystar = myfun(lambda1); No newline at end of file |
|
|
307 | 307 | Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G)); No newline at end of file |
|
|
308 | 308 | ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda); No newline at end of file |
|
|
309 | 309 | es = np.dot(Hr, f_maxent) - gnz; # should be same as ep No newline at end of file |
|
|
310 | 310 | chi2 = np.sum((es/sigma)**2); No newline at end of file |
|
|
311 | 311 | No newline at end of file |
|
|
312 | 312 | |
|
|
313 | No newline at end of file | |
|
|
313 | # CS inversion using irls ######################## No newline at end of file | |
|
|
314 | 314 | No newline at end of file |
|
|
315 | 315 | # (Use Nr, thetar, gnz, and Hr from MaxEnt above) No newline at end of file |
|
|
316 | 316 | |
|
|
317 | No newline at end of file | |
|
|
317 | Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?) No newline at end of file | |
|
|
No newline at end of file | ||
|
|
318 | ||
|
|
No newline at end of file | ||
|
|
319 | print 'FINALLY!' | |
|
|
No newline at end of file | ||
|
|
320 | print Psi.shape No newline at end of file | |
|
|
318 | 321 | No newline at end of file |
|
|
319 | 322 | # REMOVE THIS?-------------------------------- No newline at end of file |
|
|
320 | 323 | #wavelet1 = pywt.Wavelet('db4') No newline at end of file |
|
|
321 | 324 | #Phi, Psi, x = wavelet1.wavefun(level=3) No newline at end of file |
|
|
322 | 325 | # -------------------------------------------- No newline at end of file |
|
|
323 | 326 | No newline at end of file |
|
|
324 | 327 | # add "sum to 1" constraint |
|
|
328 | No newline at end of file | |
|
|
325 | H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); | |
|
|
No newline at end of file | ||
|
|
329 | No newline at end of file | |
|
|
326 | N_temp = np.array([[Nr/Nt]]); | |
|
|
No newline at end of file | ||
|
|
330 | No newline at end of file | |
|
|
327 | g2 = np.concatenate( (gnz, N_temp), axis=0 ); | |
|
|
No newline at end of file | ||
|
|
331 | No newline at end of file | |
|
|
328 | H2 = H2.T.conj(); | |
|
|
No newline at end of file | ||
|
|
332 | No newline at end of file | |
|
|
329 | ||
|
|
No newline at end of file | ||
|
|
333 | No newline at end of file | |
|
|
330 | print 'H2 shape', H2.shape | |
|
|
No newline at end of file | ||
|
|
334 | No newline at end of file | |
|
|
331 | print 'Psi shape', Psi.shape | |
|
|
No newline at end of file | ||
|
|
335 | No newline at end of file | |
|
|
332 | ||
|
|
No newline at end of file | ||
|
|
336 | No newline at end of file | |
|
|
333 | s = irls_dn2(H2*Psi,g2,0.5,G); No newline at end of file | |
|
|
334 | 337 | # f_cs = Psi*s; No newline at end of file |
|
|
335 | 338 | # No newline at end of file |
|
|
336 | 339 | # # plot No newline at end of file |
|
|
337 | 340 | # plot(thetar,f_cs,'r.-'); No newline at end of file |
|
|
338 | 341 | # hold on; No newline at end of file |
|
|
339 | 342 | # plot(thetat,fact,'k-'); No newline at end of file |
|
|
340 | 343 | # hold off; No newline at end of file |
|
|
341 | 344 | No newline at end of file |
|
|
342 | 345 | No newline at end of file |
|
|
343 | 346 | # # # Scaling and shifting No newline at end of file |
|
|
344 | 347 | # # # Only necessary for capon solution No newline at end of file |
|
|
345 | 348 | No newline at end of file |
|
|
346 | 349 | No newline at end of file |
|
|
347 | 350 | f_capon = f_capon/np.max(f_capon)*np.max(fact); No newline at end of file |
|
|
348 | 351 | No newline at end of file |
|
|
349 | 352 | No newline at end of file |
|
|
350 | 353 | ### analyze stuff ###################### No newline at end of file |
|
|
351 | 354 | # calculate MSE No newline at end of file |
|
|
352 | 355 | rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2)); No newline at end of file |
|
|
353 | 356 | rmse_capon = np.sqrt(np.mean((f_capon - factr)**2)); No newline at end of file |
|
|
354 | 357 | rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2)); No newline at end of file |
|
|
355 | 358 | #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2)); No newline at end of file |
|
|
356 | 359 | No newline at end of file |
|
|
357 | 360 | No newline at end of file |
|
|
358 | 361 | relrmse_fourier = rmse_fourier / np.linalg.norm(fact); No newline at end of file |
|
|
359 | 362 | relrmse_capon = rmse_capon / np.linalg.norm(fact); No newline at end of file |
|
|
360 | 363 | relrmse_maxent = rmse_maxent / np.linalg.norm(fact); No newline at end of file |
|
|
361 | 364 | #relrmse_cs = rmse_cs / np.norm(fact); No newline at end of file |
|
|
362 | 365 | No newline at end of file |
|
|
363 | 366 | # To be able to perform dot product (align matrices) done below within the dot calculations No newline at end of file |
|
|
364 | 367 | No newline at end of file |
|
|
365 | 368 | No newline at end of file |
|
|
366 | 369 | #f_fourier = f_fourier.T.conj() No newline at end of file |
|
|
367 | 370 | #f_capon = f_capon.T.conj() No newline at end of file |
|
|
368 | 371 | #f_maxent = f_maxent.T.conj() No newline at end of file |
|
|
369 | 372 | No newline at end of file |
|
|
370 | 373 | #factr = factr.T.conj() No newline at end of file |
|
|
371 | 374 | No newline at end of file |
|
|
372 | 375 | # calculate correlation No newline at end of file |
|
|
373 | 376 | No newline at end of file |
|
|
374 | 377 | corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr)); No newline at end of file |
|
|
375 | 378 | corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr)); No newline at end of file |
|
|
376 | 379 | corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr)); No newline at end of file |
|
|
377 | 380 | #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr)); No newline at end of file |
|
|
378 | 381 | No newline at end of file |
|
|
379 | 382 | No newline at end of file |
|
|
380 | 383 | # calculate centered correlation No newline at end of file |
|
|
381 | 384 | f0 = factr - np.mean(factr); No newline at end of file |
|
|
382 | 385 | f1 = f_fourier - np.mean(f_fourier); No newline at end of file |
|
|
383 | 386 | No newline at end of file |
|
|
384 | 387 | corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file |
|
|
385 | 388 | f1 = f_capon - np.mean(f_capon); No newline at end of file |
|
|
386 | 389 | corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file |
|
|
387 | 390 | f1 = f_maxent - np.mean(f_maxent); No newline at end of file |
|
|
388 | 391 | corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); No newline at end of file |
|
|
389 | 392 | #f1 = f_cs - mean(f_cs); No newline at end of file |
|
|
390 | 393 | #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1)); No newline at end of file |
|
|
391 | 394 | No newline at end of file |
|
|
392 | 395 | No newline at end of file |
|
|
393 | 396 | No newline at end of file |
|
|
394 | 397 | # # # plot stuff ######################### No newline at end of file |
|
|
395 | 398 | No newline at end of file |
|
|
396 | 399 | #---- Capon---- No newline at end of file |
|
|
397 | 400 | plt.figure(4) No newline at end of file |
|
|
398 | 401 | plt.subplot(2, 1, 1) No newline at end of file |
|
|
399 | 402 | plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon') No newline at end of file |
|
|
400 | 403 | plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') No newline at end of file |
|
|
401 | 404 | plt.ylabel('Power (arbitrary units)') No newline at end of file |
|
|
402 | 405 | plt.legend(loc='upper right') No newline at end of file |
|
|
403 | 406 | No newline at end of file |
|
|
404 | 407 | # formatting y-axis No newline at end of file |
|
|
405 | 408 | locs,labels = plt.yticks() No newline at end of file |
|
|
406 | 409 | plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) No newline at end of file |
|
|
407 | 410 | plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) No newline at end of file |
|
|
408 | 411 | No newline at end of file |
|
|
409 | 412 | No newline at end of file |
|
|
410 | 413 | #---- MaxEnt---- No newline at end of file |
|
|
411 | 414 | plt.subplot(2, 1, 2) No newline at end of file |
|
|
412 | 415 | plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt') No newline at end of file |
|
|
413 | 416 | plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') No newline at end of file |
|
|
414 | 417 | plt.ylabel('Power (arbitrary units)') No newline at end of file |
|
|
415 | 418 | plt.legend(loc='upper right') No newline at end of file |
|
|
416 | 419 | No newline at end of file |
|
|
417 | 420 | # formatting y-axis No newline at end of file |
|
|
418 | 421 | locs,labels = plt.yticks() No newline at end of file |
|
|
419 | 422 | plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) No newline at end of file |
|
|
420 | 423 | plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) No newline at end of file |
|
|
421 | 424 | No newline at end of file |
|
|
422 | 425 | plt.show() No newline at end of file |
|
|
423 | 426 | No newline at end of file |
|
|
424 | 427 | No newline at end of file |
|
|
425 | 428 | # # PLOT PARA COMPRESSED SENSING No newline at end of file |
|
|
426 | 429 | # # No newline at end of file |
|
|
427 | 430 | # # subplot(3,1,3); No newline at end of file |
|
|
428 | 431 | # # plot(180/pi*thetar,f_cs,'r-'); No newline at end of file |
|
|
429 | 432 | # # hold on; No newline at end of file |
|
|
430 | 433 | # # plot(180/pi*thetat,fact,'k--'); No newline at end of file |
|
|
431 | 434 | # # hold off; No newline at end of file |
|
|
432 | 435 | # # ylim([min(f_cs) 1.1*max(fact)]); No newline at end of file |
|
|
433 | 436 | # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs)); No newline at end of file |
|
|
434 | 437 | # # # title 'Compressed Sensing - Debauchies Wavelets' No newline at end of file |
|
|
435 | 438 | # # xlabel 'Degrees' No newline at end of file |
|
|
436 | 439 | # # ylabel({'Power';'(arbitrary units)'}) No newline at end of file |
|
|
437 | 440 | # # legend('Comp. Sens.','Truth'); No newline at end of file |
|
|
438 | 441 | # # No newline at end of file |
|
|
439 | 442 | # # # set(gcf,'Position',[749 143 528 881]); # CSL No newline at end of file |
|
|
440 | 443 | # # # set(gcf,'Position',[885 -21 528 673]); # macbook No newline at end of file |
|
|
441 | 444 | # # pause(0.01); No newline at end of file |
|
|
442 | 445 | No newline at end of file |
|
|
443 | 446 | No newline at end of file |
|
|
444 | 447 | # # Store Results No newline at end of file |
|
|
445 | 448 | corr[0, snri, Ni] = corr_fourier; No newline at end of file |
|
|
446 | 449 | corr[1, snri, Ni] = corr_capon; No newline at end of file |
|
|
447 | 450 | corr[2, snri, Ni] = corr_maxent; No newline at end of file |
|
|
448 | 451 | #corr[3, snri, Ni] = corr_cs; No newline at end of file |
|
|
449 | 452 | No newline at end of file |
|
|
450 | 453 | rmse[0,snri,Ni] = relrmse_fourier; No newline at end of file |
|
|
451 | 454 | rmse[1,snri,Ni] = relrmse_capon; No newline at end of file |
|
|
452 | 455 | rmse[2,snri,Ni] = relrmse_maxent; No newline at end of file |
|
|
453 | 456 | #rmse[3,snri,Ni] = relrmse_cs; No newline at end of file |
|
|
454 | 457 | No newline at end of file |
|
|
455 | 458 | corrc[0,snri,Ni] = corrc_fourier; No newline at end of file |
|
|
456 | 459 | corrc[1,snri,Ni] = corrc_capon; No newline at end of file |
|
|
457 | 460 | corrc[2,snri,Ni] = corrc_maxent; No newline at end of file |
|
|
458 | 461 | #corrc[3,snri,Ni] = corrc_cs; No newline at end of file |
|
|
459 | 462 | No newline at end of file |
|
|
460 | 463 | No newline at end of file |
|
|
461 | 464 | print 'Capon:\t', elapsed_time_capon, 'sec'; No newline at end of file |
|
|
462 | 465 | print 'Maxent:\t',elapsed_time_maxent, 'sec'; No newline at end of file |
|
|
463 | 466 | #print 'CS:\t%3.3f sec\n',elapsed_time_cs; No newline at end of file |
|
|
464 | 467 | No newline at end of file |
|
|
465 | 468 | print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN); No newline at end of file |
|
|
466 | 469 | No newline at end of file |
|
|
467 | 470 | print corr No newline at end of file |
|
|
468 | 471 | No newline at end of file |
| @@ -1,59 +1,72 | |||
|
|
1 | 1 | ''' No newline at end of file |
|
|
2 | 2 | Created on May 29, 2014 No newline at end of file |
|
|
3 | 3 | No newline at end of file |
|
|
4 | 4 | @author: Yolian Amaro No newline at end of file |
|
|
5 | 5 | ''' No newline at end of file |
|
|
6 | 6 | No newline at end of file |
|
|
7 | 7 | #from sp import multirate |
|
|
8 | No newline at end of file | |
|
|
8 | import cshift No newline at end of file | |
|
|
9 | 9 | from multirate import upfirdn No newline at end of file |
|
|
10 | 10 | No newline at end of file |
|
|
11 | 11 | def afb(x, af): No newline at end of file |
|
|
12 | 12 | No newline at end of file |
|
|
13 | 13 | # Analysis filter bank No newline at end of file |
|
|
14 | 14 | # No newline at end of file |
|
|
15 | 15 | # USAGE: No newline at end of file |
|
|
16 | 16 | # [lo, hi] = afb(x, af) No newline at end of file |
|
|
17 | 17 | # INPUT: No newline at end of file |
|
|
18 | 18 | # x - N-point vector, where No newline at end of file |
|
|
19 | 19 | # 1) N is even No newline at end of file |
|
|
20 | 20 | # 2) N >= length(af) No newline at end of file |
|
|
21 | 21 | # af - analysis filters No newline at end of file |
|
|
22 | 22 | # af(:, 1) - lowpass filter (even length) No newline at end of file |
|
|
23 | 23 | # af(:, 2) - highpass filter (even length) No newline at end of file |
|
|
24 | 24 | # OUTPUT: No newline at end of file |
|
|
25 | 25 | # lo - Low frequecy output No newline at end of file |
|
|
26 | 26 | # hi - High frequency output No newline at end of file |
|
|
27 | 27 | # EXAMPLE: No newline at end of file |
|
|
28 | 28 | # [af, sf] = farras; No newline at end of file |
|
|
29 | 29 | # x = rand(1,64); No newline at end of file |
|
|
30 | 30 | # [lo, hi] = afb(x, af); No newline at end of file |
|
|
31 | 31 | # y = sfb(lo, hi, sf); No newline at end of file |
|
|
32 | 32 | # err = x - y; No newline at end of file |
|
|
33 | 33 | # max(abs(err)) No newline at end of file |
|
|
34 | 34 | # No newline at end of file |
|
|
35 | 35 | # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY No newline at end of file |
|
|
36 | 36 | # http://taco.poly.edu/WaveletSoftware/ No newline at end of file |
|
|
37 | 37 | No newline at end of file |
|
|
38 | 38 | N = x.size; |
|
|
39 | No newline at end of file | |
|
|
39 | L = (af).size/2; | |
|
|
No newline at end of file | ||
|
|
40 | No newline at end of file | |
|
|
40 | x = cshift(x,-L); No newline at end of file | |
|
|
No newline at end of file | ||
|
|
41 | #print 'L', L | |
|
|
No newline at end of file | ||
|
|
42 | x = cshift(x,-(L-1)); | |
|
|
No newline at end of file | ||
|
|
43 | ||
|
|
No newline at end of file | ||
|
|
44 | # print 'afb x', x.shape | |
|
|
No newline at end of file | ||
|
|
45 | # print 'af[:,0]',af[:,0].shape | |
|
|
No newline at end of file | ||
|
|
46 | # print 'af[:,1]',af[:,1].shape | |
|
|
No newline at end of file | ||
|
|
47 | # print '-----------------------' No newline at end of file | |
|
|
41 | 48 | No newline at end of file |
|
|
42 | 49 | # lowpass filter No newline at end of file |
|
|
43 | 50 | lo = upfirdn(x, af[:,0], 1, 2); No newline at end of file |
|
|
44 | 51 | No newline at end of file |
|
|
45 | 52 | # VERIFY THIS!!!!!!!!!!!! No newline at end of file |
|
|
46 | 53 | for i in range(0, L): |
|
|
54 | No newline at end of file | |
|
|
47 | lo[i] = lo[N/2+[i]] + lo[i]; No newline at end of file | |
|
|
48 | 55 | |
|
|
56 | No newline at end of file | |
|
|
49 | lo = lo[1:N/2]; No newline at end of file | |
|
|
No newline at end of file | ||
|
|
57 | lo = lo[0:N/2]; No newline at end of file | |
|
|
50 | 58 | No newline at end of file |
|
|
51 | 59 | # highpass filter |
|
|
60 | No newline at end of file | |
|
|
52 | hi = upfirdn(x, af[:,2], 1, 2); No newline at end of file | |
|
|
53 | 61 | No newline at end of file |
|
|
54 | 62 | for j in range(0, L): |
|
|
63 | No newline at end of file | |
|
|
55 | hi[j] = hi(N/2+[j]) + hi[j]; No newline at end of file | |
|
|
56 | 64 | |
|
|
65 | No newline at end of file | |
|
|
57 | hi = hi[1:N/2]; No newline at end of file | |
|
|
No newline at end of file | ||
|
|
66 | hi = hi[0:N/2]; | |
|
|
No newline at end of file | ||
|
|
67 | ||
|
|
No newline at end of file | ||
|
|
68 | # Reshape from 1D to 2D | |
|
|
No newline at end of file | ||
|
|
69 | lo = lo.reshape(1, lo.size) No newline at end of file | |
|
|
58 | 70 | No newline at end of file |
|
|
59 | 71 | return lo, hi No newline at end of file |
|
|
60 | 72 | |
| @@ -1,29 +1,34 | |||
|
|
1 | 1 | ''' No newline at end of file |
|
|
2 | 2 | Created on May 30, 2014 No newline at end of file |
|
|
3 | 3 | No newline at end of file |
|
|
4 | 4 | @author: Yolian Amaro No newline at end of file |
|
|
5 | 5 | ''' No newline at end of file |
|
|
6 | 6 | No newline at end of file |
|
|
7 | 7 | import numpy as np No newline at end of file |
|
|
8 | 8 | No newline at end of file |
|
|
9 | 9 | def cshift(x, m): No newline at end of file |
|
|
10 | 10 | No newline at end of file |
|
|
11 | 11 | # Circular Shift No newline at end of file |
|
|
12 | 12 | # No newline at end of file |
|
|
13 | 13 | # USAGE: No newline at end of file |
|
|
14 | 14 | # y = cshift(x, m) No newline at end of file |
|
|
15 | 15 | # INPUT: No newline at end of file |
|
|
16 | 16 | # x - N-point vector |
|
|
17 | No newline at end of file | |
|
|
17 | # m - amount of shift No newline at end of file | |
|
|
18 | 18 | # OUTPUT: No newline at end of file |
|
|
19 | 19 | # y - vector x will be shifted by m samples to the left No newline at end of file |
|
|
20 | 20 | # No newline at end of file |
|
|
21 | 21 | # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY No newline at end of file |
|
|
22 | 22 | # http://taco.poly.edu/WaveletSoftware/ No newline at end of file |
|
|
23 | 23 | No newline at end of file |
|
|
24 | 24 | N = x.size; |
|
|
25 | No newline at end of file | |
|
|
25 | n = np.arange(N-1); No newline at end of file | |
|
|
26 | 26 |
n = np. |
|
|
27 | No newline at end of file | |
|
|
27 | y = x[n]; No newline at end of file | |
|
|
No newline at end of file | ||
|
|
28 | ||
|
|
No newline at end of file | ||
|
|
29 | print x.shape | |
|
|
No newline at end of file | ||
|
|
30 | ||
|
|
No newline at end of file | ||
|
|
31 | y = x[0,n]; No newline at end of file | |
|
|
28 | 32 | No newline at end of file |
|
|
29 | 33 | return y No newline at end of file |
|
|
30 | 34 | L30: rhodecode diff rendering error |
| @@ -1,38 +1,41 | |||
|
|
1 | 1 | ''' No newline at end of file |
|
|
2 | 2 | Created on May 26, 2014 No newline at end of file |
|
|
3 | 3 | No newline at end of file |
|
|
4 | 4 | @author: Yolian Amaro No newline at end of file |
|
|
5 | 5 | ''' No newline at end of file |
|
|
6 | 6 | No newline at end of file |
|
|
7 | 7 | import numpy as np |
|
|
8 | No newline at end of file | |
|
|
8 | import FSfarras | |
|
|
No newline at end of file | ||
|
|
9 | No newline at end of file | |
|
|
9 | import dualfilt1 No newline at end of file | |
|
|
No newline at end of file | ||
|
|
10 | from dualtree import * | |
|
|
No newline at end of file | ||
|
|
11 | from idualtree import * No newline at end of file | |
|
|
10 | 12 | No newline at end of file |
|
|
11 | 13 | def deb4_basis(N): No newline at end of file |
|
|
12 | 14 | No newline at end of file |
|
|
13 | 15 | Psi = np.zeros(shape=(N,2*N+1)); No newline at end of file |
|
|
14 | 16 | idx = 1; No newline at end of file |
|
|
15 | 17 | No newline at end of file |
|
|
16 | 18 | J = 4; |
|
|
19 | No newline at end of file | |
|
|
17 | [Faf, Fsf] = FSfarras; | |
|
|
No newline at end of file | ||
|
|
20 | No newline at end of file | |
|
|
18 | [af, sf] = dualfilt1; No newline at end of file | |
|
|
19 | 21 | No newline at end of file |
|
|
20 | 22 | # compute transform of zero vector No newline at end of file |
|
|
21 | 23 | x = np.zeros(shape=(1,N)); |
|
|
24 | No newline at end of file | |
|
|
22 | #w = dualtree(x, J, Faf, af); | |
|
|
No newline at end of file | ||
|
|
25 | No newline at end of file | |
|
|
23 | # # | |
|
|
No newline at end of file | ||
|
|
26 | No newline at end of file | |
|
|
24 | # # # Uses both real and imaginary wavelets | |
|
|
No newline at end of file | ||
|
|
27 | No newline at end of file | |
|
|
25 | # # for i in range (1, J+1): | |
|
|
No newline at end of file | ||
|
|
28 | No newline at end of file | |
|
|
26 | # # for j in range (1, 2): | |
|
|
No newline at end of file | ||
|
|
29 | No newline at end of file | |
|
|
27 | # # for k in range (1, (w[i][j]).size): | |
|
|
No newline at end of file | ||
|
|
30 | No newline at end of file | |
|
|
28 | # # w[i][j](k) = 1; | |
|
|
No newline at end of file | ||
|
|
31 | No newline at end of file | |
|
|
29 | # # y = idualtree(w, J, Fsf, sf); | |
|
|
No newline at end of file | ||
|
|
32 | No newline at end of file | |
|
|
30 | # # w[i][j](k) = 0; | |
|
|
No newline at end of file | ||
|
|
33 | No newline at end of file | |
|
|
31 | # # # store it | |
|
|
No newline at end of file | ||
|
|
34 | No newline at end of file | |
|
|
32 | # # Psi(:,idx) = y.T.conj(); | |
|
|
No newline at end of file | ||
|
|
35 | No newline at end of file | |
|
|
33 | # # idx = idx + 1; | |
|
|
No newline at end of file | ||
|
|
36 | No newline at end of file | |
|
|
34 | # # | |
|
|
No newline at end of file | ||
|
|
37 | No newline at end of file | |
|
|
35 | # # # Add uniform vector (seems to be useful if there's a background | |
|
|
No newline at end of file | ||
|
|
38 | No newline at end of file | |
|
|
36 | # # Psi(:,2*N+1) = 1/np.sqrt(N); | |
|
|
No newline at end of file | ||
|
|
39 | No newline at end of file | |
|
|
37 | # | |
|
|
No newline at end of file | ||
|
|
40 | No newline at end of file | |
|
|
38 | # return Psi No newline at end of file | |
|
|
No newline at end of file | ||
|
|
41 | return Psi No newline at end of file | |
| @@ -1,63 +1,76 | |||
|
|
1 | 1 | ''' No newline at end of file |
|
|
2 | 2 | Created on May 29, 2014 No newline at end of file |
|
|
3 | 3 | No newline at end of file |
|
|
4 | 4 | @author: Yolian Amaro No newline at end of file |
|
|
5 | 5 | ''' No newline at end of file |
|
|
6 | 6 | No newline at end of file |
|
|
7 | 7 | import numpy as np |
|
|
8 | No newline at end of file | |
|
|
8 | import afb No newline at end of file | |
|
|
9 | 9 | No newline at end of file |
|
|
10 | 10 | def dualtree(x, J, Faf, af): No newline at end of file |
|
|
11 | 11 | No newline at end of file |
|
|
12 | 12 | # Dual-tree Complex Discrete Wavelet Transform No newline at end of file |
|
|
13 | 13 | # No newline at end of file |
|
|
14 | 14 | # USAGE: No newline at end of file |
|
|
15 | 15 | # w = dualtree(x, J, Faf, af) No newline at end of file |
|
|
16 | 16 | # INPUT: No newline at end of file |
|
|
17 | 17 | # x - N-point vector No newline at end of file |
|
|
18 | 18 | # 1) N is divisible by 2^J No newline at end of file |
|
|
19 | 19 | # 2) N >= 2^(J-1)*length(af) No newline at end of file |
|
|
20 | 20 | # J - number of stages No newline at end of file |
|
|
21 | 21 | # Faf - filters for the first stage No newline at end of file |
|
|
22 | 22 | # af - filters for the remaining stages No newline at end of file |
|
|
23 | 23 | # OUTPUT: No newline at end of file |
|
|
24 | 24 | # w - DWT coefficients No newline at end of file |
|
|
25 | 25 | # w{j}{1}, j = 1..J - real part No newline at end of file |
|
|
26 | 26 | # w{j}{2}, j = 1..J - imaginary part No newline at end of file |
|
|
27 | 27 | # w{J+1}{d} - lowpass coefficients, d = 1,2 No newline at end of file |
|
|
28 | 28 | # EXAMPLE: No newline at end of file |
|
|
29 | 29 | # x = rand(1, 512); No newline at end of file |
|
|
30 | 30 | # J = 4; No newline at end of file |
|
|
31 | 31 | # [Faf, Fsf] = FSfarras; No newline at end of file |
|
|
32 | 32 | # [af, sf] = dualfilt1; No newline at end of file |
|
|
33 | 33 | # w = dualtree(x, J, Faf, af); No newline at end of file |
|
|
34 | 34 | # y = idualtree(w, J, Fsf, sf); No newline at end of file |
|
|
35 | 35 | # err = x - y; No newline at end of file |
|
|
36 | 36 | # max(abs(err)) No newline at end of file |
|
|
37 | 37 | # No newline at end of file |
|
|
38 | 38 | # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY No newline at end of file |
|
|
39 | 39 | # http://taco.poly.edu/WaveletSoftware/ No newline at end of file |
|
|
40 | 40 | No newline at end of file |
|
|
41 | # ---------Trees Structure---------------# | |
|
|
No newline at end of file | ||
|
|
42 | # w [ 0 1 2 .... J ] # | |
|
|
No newline at end of file | ||
|
|
43 | # | | | | # | |
|
|
No newline at end of file | ||
|
|
44 | # [0 1] [0 1] [0 1] [0 1] # | |
|
|
No newline at end of file | ||
|
|
45 | #----------------------------------------# | |
|
|
No newline at end of file | ||
|
|
46 | No newline at end of file | |
|
|
41 | 47 | # normalization No newline at end of file |
|
|
42 | 48 | x = x/np.sqrt(2); No newline at end of file |
|
|
43 | 49 | |
|
|
50 | No newline at end of file | |
|
|
44 | w = np.zeros(shape=(J,2)) ### VERIFY THIS DEFINITION No newline at end of file | |
|
|
No newline at end of file | ||
|
|
51 | w = np.zeros(shape=(J+1), dtype=object) | |
|
|
No newline at end of file | ||
|
|
52 | ||
|
|
No newline at end of file | ||
|
|
53 | for j in range (0, w.size): | |
|
|
No newline at end of file | ||
|
|
54 | w[j] = np.zeros(shape=(J+1), dtype=object) No newline at end of file | |
|
|
45 | 55 | No newline at end of file |
|
|
46 | 56 | # Tree 1 |
|
|
57 | No newline at end of file | |
|
|
47 | [x1,w[1,0]] = afb(x, Faf[0,1]); # w{1}{1} No newline at end of file | |
|
|
48 | 58 | |
|
|
No newline at end of file | ||
|
|
49 | for j in range (2,J): | |
|
|
No newline at end of file | ||
|
|
50 | [x1,w[j,0]] = afb(x1, af[0,1]); #check this No newline at end of file | |
|
|
51 | 59 |
|
|
|
60 | No newline at end of file | |
|
|
52 | w[J+1,1] = x1; No newline at end of file | |
|
|
No newline at end of file | ||
|
|
61 | [x1,w[j][0]] = afb(x1, af[0,0]); ### or 0,1???? | |
|
|
No newline at end of file | ||
|
|
62 | ||
|
|
No newline at end of file | ||
|
|
63 | ||
|
|
No newline at end of file | ||
|
|
64 | ||
|
|
No newline at end of file | ||
|
|
65 | w[J][0] = x1; No newline at end of file | |
|
|
53 | 66 | No newline at end of file |
|
|
54 | 67 | # Tree 2 |
|
|
68 | No newline at end of file | |
|
|
55 | [x2,w[1,2]] = afb(x, Faf[0,1]); No newline at end of file | |
|
|
56 | 69 | |
|
|
70 | No newline at end of file | |
|
|
57 | for j in range (2,J): | |
|
|
No newline at end of file | ||
|
|
71 | No newline at end of file | |
|
|
58 | [x2,w[j,2]] = afb(x2, af[0,1]); No newline at end of file | |
|
|
59 | 72 | |
|
|
73 | No newline at end of file | |
|
|
60 | w[J+1,2] = x2; No newline at end of file | |
|
|
61 | 74 | No newline at end of file |
|
|
62 | 75 | return w No newline at end of file |
|
|
63 | 76 | No newline at end of file |
General Comments 0
You need to be logged in to leave comments.
Login now