##// 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 26 # Translated to Python by Yolian Amaro
27 27
28 28
29
30 29 a1 = np.array( [
31 30 [ 0, 0],
32 31 [-0.08838834764832, -0.01122679215254],
@@ -53,7 +52,6
53 52 [ 0, -0.01122679215254]
54 53 ]);
55 54
56 #print a2.shape
57 55
58 56 af = np.array([ [a1,a2] ], dtype=object)
59 57
Show file before | Show file after | Hide comments
@@ -310,11 +310,14
310 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 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 322 # REMOVE THIS?--------------------------------
320 323 #wavelet1 = pywt.Wavelet('db4')
@@ -322,15 +325,15
322 325 # --------------------------------------------
323 326
324 327 # add "sum to 1" constraint
325 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
326 N_temp = np.array([[Nr/Nt]]);
327 g2 = np.concatenate( (gnz, N_temp), axis=0 );
328 H2 = H2.T.conj();
329
330 print 'H2 shape', H2.shape
331 print 'Psi shape', Psi.shape
332
333 s = irls_dn2(H2*Psi,g2,0.5,G);
328 # H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
329 # N_temp = np.array([[Nr/Nt]]);
330 # g2 = np.concatenate( (gnz, N_temp), axis=0 );
331 # H2 = H2.T.conj();
332 #
333 # print 'H2 shape', H2.shape
334 # print 'Psi shape', Psi.shape
335 #
336 # s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
334 337 # f_cs = Psi*s;
335 338 #
336 339 # # plot
Show file before | Show file after | Hide comments
@@ -5,7 +5,7
5 5 '''
6 6
7 7 #from sp import multirate
8 import cshift
8 from cshift import *
9 9 from multirate import upfirdn
10 10
11 11 def afb(x, af):
@@ -36,24 +36,37
36 36 # http://taco.poly.edu/WaveletSoftware/
37 37
38 38 N = x.size;
39 L = (af).size/2;
40 x = cshift(x,-L);
39 L = (af).size/4; #L should be = 5
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 49 # lowpass filter
43 50 lo = upfirdn(x, af[:,0], 1, 2);
44 51
52
45 53 # VERIFY THIS!!!!!!!!!!!!
46 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 60 # highpass filter
52 hi = upfirdn(x, af[:,2], 1, 2);
61 hi = upfirdn(x, af[:,1], 1, 2);
53 62
54 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 72 return lo, hi
Show file before | Show file after | Hide comments
@@ -14,16 +14,21
14 14 # y = cshift(x, m)
15 15 # INPUT:
16 16 # x - N-point vector
17 # m - amount of shift
17 # m - amount of shift (pos=left, neg=right)
18 18 # OUTPUT:
19 19 # y - vector x will be shifted by m samples to the left
20 20 #
21 21 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
22 22 # http://taco.poly.edu/WaveletSoftware/
23
23 24
24 25 N = x.size;
25 n = np.arange(N-1);
26 n = np.arange(N);
26 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 34 return y No newline at end of file
Show file before | Show file after | Hide comments
@@ -5,8 +5,10
5 5 '''
6 6
7 7 import numpy as np
8 import FSfarras
9 import dualfilt1
8 from FSfarras import *
9 from dualfilt1 import *
10 from dualtree import *
11 from idualtree import *
10 12
11 13 def deb4_basis(N):
12 14
@@ -14,25 +16,26
14 16 idx = 1;
15 17
16 18 J = 4;
17 [Faf, Fsf] = FSfarras;
18 [af, sf] = dualfilt1;
19 [Faf, Fsf] = FSfarras();
20 [af, sf] = dualfilt1();
19 21
20 22 # compute transform of zero vector
21 23 x = np.zeros(shape=(1,N));
22 #w = dualtree(x, J, Faf, af);
23 # #
24 # # # Uses both real and imaginary wavelets
25 # # for i in range (1, J+1):
26 # # for j in range (1, 2):
27 # # for k in range (1, (w[i][j]).size):
28 # # w[i][j](k) = 1;
29 # # y = idualtree(w, J, Fsf, sf);
30 # # w[i][j](k) = 0;
31 # # # store it
32 # # Psi(:,idx) = y.T.conj();
33 # # idx = idx + 1;
34 # #
35 # # # Add uniform vector (seems to be useful if there's a background
36 # # Psi(:,2*N+1) = 1/np.sqrt(N);
37 #
38 # return Psi No newline at end of file
24 w = dualtree(x, J, Faf, af);
25
26
27 # Uses both real and imaginary wavelets
28 for i in range (0, J):
29 for j in range (0, 1):
30 for k in range (0, (w[i][j]).size):
31 w[i][j][0,k] = 1;
32 y = idualtree(w, J, Fsf, sf);
33 w[i][j][0,k] = 0;
34 # store it
35 Psi[:,idx] = y.T.conj();
36 idx = idx + 1;
37
38 # Add uniform vector (seems to be useful if there's a background
39 Psi[:,2*N+1] = 1/np.sqrt(N);
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 7 import numpy as np
8 import afb
8 from afb import *
9 9
10 10 def dualtree(x, J, Faf, af):
11 11
@@ -38,26 +38,39
38 38 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
39 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 47 # normalization
42 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 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 67 # Tree 2
55 [x2,w[1,2]] = afb(x, Faf[0,1]);
56
57 for j in range (2,J):
58 [x2,w[j,2]] = afb(x2, af[0,1]);
59
60 w[J+1,2] = x2;
68 [x2,w[0][1]] = afb(x, Faf[0,1]);
69
70 for j in range (1,J):
71 [x2,w[j][1]] = afb(x2, af[0,1]);
72
73 w[J][1] = x2;
61 74
62 75 return w
63 76
General Comments 0
You need to be logged in to leave comments. Login now