##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
Fixed errors, implemented sfb and idualtree classes
yamaro -
r19:20
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,46
1 '''
2 Created on Jun 5, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 from sfb import *
8
9 def idualtree(w, J, Fsf, sf):
10
11 # Inverse Dual-tree Complex DWT
12 #
13 # USAGE:
14 # y = idualtree(w, J, Fsf, sf)
15 # INPUT:
16 # w - DWT coefficients
17 # J - number of stages
18 # Fsf - synthesis filters for the last stage
19 # sf - synthesis filters for preceeding stages
20 # OUTUT:
21 # y - output signal
22 # See dualtree
23 #
24 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
25 # http://taco.poly.edu/WaveletSoftware/
26
27 # Tree 1
28 y1 = w[J][0];
29
30 for j in range (J-1, 0, -1):
31 y1 = sfb(y1, w[j][0], sf[0,0]);
32
33 y1 = sfb(y1, w[0][0], Fsf[0,0]);
34
35 # Tree 2
36 y2 = w[J][1];
37
38 for j in range (J-1, 0, -1):
39 y2 = sfb(y2, w[j][2], sf[0,1]);
40
41 y2 = sfb(y2, w[0][1], Fsf[0,1]);
42
43 # normalization
44 y = (y1 + y2)/np.sqrt(2);
45
46 return y
Show file before | Show file after | Hide comments
@@ -0,0 +1,68
1 '''
2 Created on Jun 5, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 from multirate import *
8 import numpy as np
9 from cshift import *
10
11 def sfb(lo, hi, sf):
12
13 # Synthesis filter bank
14 #
15 # USAGE:
16 # y = sfb(lo, hi, sf)
17 # INPUT:
18 # lo - low frqeuency input
19 # hi - high frequency input
20 # sf - synthesis filters
21 # sf(:, 1) - lowpass filter (even length)
22 # sf(:, 2) - highpass filter (even length)
23 # OUTPUT:
24 # y - output signal
25 # See also afb
26 #
27 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
28 # http://taco.poly.edu/WaveletSoftware/
29
30 N = 2*lo.size;
31 L = sf.size/2;
32 #print 'N', N
33 #print 'sf', sf
34
35
36 #print 'sf[:,0]', sf[:,0].shape
37 #print 'sf[:,1]', sf[:,1].shape
38 #print 'sbf hi', hi.shape
39
40
41
42 # Need to change format for upfirdn funct:
43 lo = lo.T.conj()
44 lo = lo.reshape(lo.size)
45
46 print 'sfb hi', hi
47
48 # Need to change format for upfirdn funct:
49 hi = hi.T.conj()
50 hi = hi.reshape(hi.size)
51
52 #hi = hi.reshape(1, hi.size)
53
54 lo = upfirdn(lo, sf[:,0], 2, 1);
55 hi = upfirdn(hi, sf[:,1], 2, 1);
56 y = lo + hi;
57 y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT
58 y = y[0:N];
59
60 print 'y en sbf\n', y.shape
61
62 y = y.reshape(1, y.size)
63 print 'y en sbf\n', y.shape
64
65 y = cshift(y, 1-L/2);
66
67 return y;
68
Show file before | Show file after | Hide comments
@@ -26,7 +26,6
26 # Translated to Python by Yolian Amaro
26 # Translated to Python by Yolian Amaro
27
27
28
28
29
30 a1 = np.array( [
29 a1 = np.array( [
31 [ 0, 0],
30 [ 0, 0],
32 [-0.08838834764832, -0.01122679215254],
31 [-0.08838834764832, -0.01122679215254],
@@ -53,7 +52,6
53 [ 0, -0.01122679215254]
52 [ 0, -0.01122679215254]
54 ]);
53 ]);
55
54
56 #print a2.shape
57
55
58 af = np.array([ [a1,a2] ], dtype=object)
56 af = np.array([ [a1,a2] ], dtype=object)
59
57
Show file before | Show file after | Hide comments
@@ -310,11 +310,14
310 chi2 = np.sum((es/sigma)**2);
310 chi2 = np.sum((es/sigma)**2);
311
311
312
312
313 # CS inversion using irls ########################
313 # CS inversion using Iteratively Reweighted Least Squares (IRLS)-------------
314
314
315 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
315 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
316
316
317 Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
317 Psi = deb4_basis(Nr); ###### REPLACED BY LINEs BELOW (?)
318
319 print 'FINALLY!'
320 print Psi.shape
318
321
319 # REMOVE THIS?--------------------------------
322 # REMOVE THIS?--------------------------------
320 #wavelet1 = pywt.Wavelet('db4')
323 #wavelet1 = pywt.Wavelet('db4')
@@ -322,15 +325,15
322 # --------------------------------------------
325 # --------------------------------------------
323
326
324 # add "sum to 1" constraint
327 # add "sum to 1" constraint
325 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
328 # H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
326 N_temp = np.array([[Nr/Nt]]);
329 # N_temp = np.array([[Nr/Nt]]);
327 g2 = np.concatenate( (gnz, N_temp), axis=0 );
330 # g2 = np.concatenate( (gnz, N_temp), axis=0 );
328 H2 = H2.T.conj();
331 # H2 = H2.T.conj();
329
332 #
330 print 'H2 shape', H2.shape
333 # print 'H2 shape', H2.shape
331 print 'Psi shape', Psi.shape
334 # print 'Psi shape', Psi.shape
332
335 #
333 s = irls_dn2(H2*Psi,g2,0.5,G);
336 # s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
334 # f_cs = Psi*s;
337 # f_cs = Psi*s;
335 #
338 #
336 # # plot
339 # # plot
Show file before | Show file after | Hide comments
@@ -5,7 +5,7
5 '''
5 '''
6
6
7 #from sp import multirate
7 #from sp import multirate
8 import cshift
8 from cshift import *
9 from multirate import upfirdn
9 from multirate import upfirdn
10
10
11 def afb(x, af):
11 def afb(x, af):
@@ -36,24 +36,37
36 # http://taco.poly.edu/WaveletSoftware/
36 # http://taco.poly.edu/WaveletSoftware/
37
37
38 N = x.size;
38 N = x.size;
39 L = (af).size/2;
39 L = (af).size/4; #L should be = 5
40 x = cshift(x,-L);
40 #print af
41 #print 'L', L
42 x = cshift(x,-(L-1));
43
44 # print 'afb x', x.shape
45 # print 'af[:,0]',af[:,0].shape
46 # print 'af[:,1]',af[:,1].shape
47 # print '-----------------------'
41
48
42 # lowpass filter
49 # lowpass filter
43 lo = upfirdn(x, af[:,0], 1, 2);
50 lo = upfirdn(x, af[:,0], 1, 2);
44
51
52
45 # VERIFY THIS!!!!!!!!!!!!
53 # VERIFY THIS!!!!!!!!!!!!
46 for i in range(0, L):
54 for i in range(0, L):
47 lo[i] = lo[N/2+[i]] + lo[i];
55 lo[i] = lo[N/2+i] + lo[i];
48
56
49 lo = lo[1:N/2];
57 lo = lo[0:N/2];
58
50
59
51 # highpass filter
60 # highpass filter
52 hi = upfirdn(x, af[:,2], 1, 2);
61 hi = upfirdn(x, af[:,1], 1, 2);
53
62
54 for j in range(0, L):
63 for j in range(0, L):
55 hi[j] = hi(N/2+[j]) + hi[j];
64 hi[j] = hi[N/2+j] + hi[j];
56
65
57 hi = hi[1:N/2];
66 hi = hi[0:N/2];
67
68 # Reshape from 1D to 2D
69 lo = lo.reshape(1, lo.size)
70 hi = hi.reshape(1, hi.size)
58
71
59 return lo, hi
72 return lo, hi
Show file before | Show file after | Hide comments
@@ -14,16 +14,21
14 # y = cshift(x, m)
14 # y = cshift(x, m)
15 # INPUT:
15 # INPUT:
16 # x - N-point vector
16 # x - N-point vector
17 # m - amount of shift
17 # m - amount of shift (pos=left, neg=right)
18 # OUTPUT:
18 # OUTPUT:
19 # y - vector x will be shifted by m samples to the left
19 # y - vector x will be shifted by m samples to the left
20 #
20 #
21 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
21 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
22 # http://taco.poly.edu/WaveletSoftware/
22 # http://taco.poly.edu/WaveletSoftware/
23
23
24
24 N = x.size;
25 N = x.size;
25 n = np.arange(N-1);
26 n = np.arange(N);
26 n = np.mod(n-m, N);
27 n = np.mod(n-m, N);
27 y = x[n];
28
29 print x.shape
30
31 y = x[0,n];
32
28
33
29 return y No newline at end of file
34 return y
Show file before | Show file after | Hide comments
@@ -5,8 +5,10
5 '''
5 '''
6
6
7 import numpy as np
7 import numpy as np
8 import FSfarras
8 from FSfarras import *
9 import dualfilt1
9 from dualfilt1 import *
10 from dualtree import *
11 from idualtree import *
10
12
11 def deb4_basis(N):
13 def deb4_basis(N):
12
14
@@ -14,25 +16,26
14 idx = 1;
16 idx = 1;
15
17
16 J = 4;
18 J = 4;
17 [Faf, Fsf] = FSfarras;
19 [Faf, Fsf] = FSfarras();
18 [af, sf] = dualfilt1;
20 [af, sf] = dualfilt1();
19
21
20 # compute transform of zero vector
22 # compute transform of zero vector
21 x = np.zeros(shape=(1,N));
23 x = np.zeros(shape=(1,N));
22 #w = dualtree(x, J, Faf, af);
24 w = dualtree(x, J, Faf, af);
23 # #
25
24 # # # Uses both real and imaginary wavelets
26
25 # # for i in range (1, J+1):
27 # Uses both real and imaginary wavelets
26 # # for j in range (1, 2):
28 for i in range (0, J):
27 # # for k in range (1, (w[i][j]).size):
29 for j in range (0, 1):
28 # # w[i][j](k) = 1;
30 for k in range (0, (w[i][j]).size):
29 # # y = idualtree(w, J, Fsf, sf);
31 w[i][j][0,k] = 1;
30 # # w[i][j](k) = 0;
32 y = idualtree(w, J, Fsf, sf);
31 # # # store it
33 w[i][j][0,k] = 0;
32 # # Psi(:,idx) = y.T.conj();
34 # store it
33 # # idx = idx + 1;
35 Psi[:,idx] = y.T.conj();
34 # #
36 idx = idx + 1;
35 # # # Add uniform vector (seems to be useful if there's a background
37
36 # # Psi(:,2*N+1) = 1/np.sqrt(N);
38 # Add uniform vector (seems to be useful if there's a background
37 #
39 Psi[:,2*N+1] = 1/np.sqrt(N);
38 # return Psi No newline at end of file
40
41 return Psi No newline at end of file
Show file before | Show file after | Hide comments
@@ -5,7 +5,7
5 '''
5 '''
6
6
7 import numpy as np
7 import numpy as np
8 import afb
8 from afb import *
9
9
10 def dualtree(x, J, Faf, af):
10 def dualtree(x, J, Faf, af):
11
11
@@ -38,26 +38,39
38 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
38 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
39 # http://taco.poly.edu/WaveletSoftware/
39 # http://taco.poly.edu/WaveletSoftware/
40
40
41 # ---------Trees Structure---------------#
42 # w [ 0 1 2 .... J ] #
43 # | | | | #
44 # [0 1] [0 1] [0 1] [0 1] #
45 #----------------------------------------#
46
41 # normalization
47 # normalization
42 x = x/np.sqrt(2);
48 x = x/np.sqrt(2);
49
43
50
44 w = np.zeros(shape=(J,2)) ### VERIFY THIS DEFINITION
51 w = np.zeros(shape=(J+1), dtype=object)
45
52
53 for j in range (0, w.size):
54 w[j] = np.zeros(shape=(J+1), dtype=object)
55
46 # Tree 1
56 # Tree 1
47 [x1,w[1,0]] = afb(x, Faf[0,1]); # w{1}{1}
57 [x1, w[0][0]] = afb(x, Faf[0,0]); # w{1}{1}
48
58
49 for j in range (2,J):
50 [x1,w[j,0]] = afb(x1, af[0,1]); #check this
51
52 w[J+1,1] = x1;
53
59
60 for j in range (1,J):
61 [x1,w[j][0]] = afb(x1, af[0,0]); ### or 0,1????
62
63
64
65 w[J][0] = x1;
66
54 # Tree 2
67 # Tree 2
55 [x2,w[1,2]] = afb(x, Faf[0,1]);
68 [x2,w[0][1]] = afb(x, Faf[0,1]);
56
69
57 for j in range (2,J):
70 for j in range (1,J):
58 [x2,w[j,2]] = afb(x2, af[0,1]);
71 [x2,w[j][1]] = afb(x2, af[0,1]);
59
72
60 w[J+1,2] = x2;
73 w[J][1] = x2;
61
74
62 return w
75 return w
63
76
General Comments 0
You need to be logged in to leave comments. Login now