Index: trunk/src/src/FSfarras.py =================================================================== diff --git a/trunk/src/src/FSfarras.py b/trunk/src/src/FSfarras.py --- a/trunk/src/src/FSfarras.py (revision 19) +++ b/trunk/src/src/FSfarras.py (revision 20) @@ -26,7 +26,6 @@ # Translated to Python by Yolian Amaro - a1 = np.array( [ [ 0, 0], [-0.08838834764832, -0.01122679215254], @@ -53,7 +52,6 @@ [ 0, -0.01122679215254] ]); - #print a2.shape af = np.array([ [a1,a2] ], dtype=object) Index: trunk/src/src/RICS_paper_compare_noise.py =================================================================== diff --git a/trunk/src/src/RICS_paper_compare_noise.py b/trunk/src/src/RICS_paper_compare_noise.py --- a/trunk/src/src/RICS_paper_compare_noise.py (revision 19) +++ b/trunk/src/src/RICS_paper_compare_noise.py (revision 20) @@ -310,11 +310,14 @@ chi2 = np.sum((es/sigma)**2); - # CS inversion using irls ######################## + # CS inversion using Iteratively Reweighted Least Squares (IRLS)------------- # (Use Nr, thetar, gnz, and Hr from MaxEnt above) - Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?) + Psi = deb4_basis(Nr); ###### REPLACED BY LINEs BELOW (?) + + print 'FINALLY!' + print Psi.shape # REMOVE THIS?-------------------------------- #wavelet1 = pywt.Wavelet('db4') @@ -322,15 +325,15 @@ # -------------------------------------------- # add "sum to 1" constraint - H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); - N_temp = np.array([[Nr/Nt]]); - g2 = np.concatenate( (gnz, N_temp), axis=0 ); - H2 = H2.T.conj(); - - print 'H2 shape', H2.shape - print 'Psi shape', Psi.shape - - s = irls_dn2(H2*Psi,g2,0.5,G); + # H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); + # N_temp = np.array([[Nr/Nt]]); + # g2 = np.concatenate( (gnz, N_temp), axis=0 ); + # H2 = H2.T.conj(); + # + # print 'H2 shape', H2.shape + # print 'Psi shape', Psi.shape + # + # s = irls_dn2(np.dot(H2,Psi),g2,0.5,G); # f_cs = Psi*s; # # # plot Index: trunk/src/src/afb.py =================================================================== diff --git a/trunk/src/src/afb.py b/trunk/src/src/afb.py --- a/trunk/src/src/afb.py (revision 19) +++ b/trunk/src/src/afb.py (revision 20) @@ -5,7 +5,7 @@ ''' #from sp import multirate -import cshift +from cshift import * from multirate import upfirdn def afb(x, af): @@ -36,24 +36,37 @@ # http://taco.poly.edu/WaveletSoftware/ N = x.size; - L = (af).size/2; - x = cshift(x,-L); + L = (af).size/4; #L should be = 5 + #print af + #print 'L', L + x = cshift(x,-(L-1)); + + # print 'afb x', x.shape + # print 'af[:,0]',af[:,0].shape + # print 'af[:,1]',af[:,1].shape + # print '-----------------------' # lowpass filter lo = upfirdn(x, af[:,0], 1, 2); + # VERIFY THIS!!!!!!!!!!!! for i in range(0, L): - lo[i] = lo[N/2+[i]] + lo[i]; + lo[i] = lo[N/2+i] + lo[i]; - lo = lo[1:N/2]; + lo = lo[0:N/2]; + # highpass filter - hi = upfirdn(x, af[:,2], 1, 2); + hi = upfirdn(x, af[:,1], 1, 2); for j in range(0, L): - hi[j] = hi(N/2+[j]) + hi[j]; + hi[j] = hi[N/2+j] + hi[j]; - hi = hi[1:N/2]; + hi = hi[0:N/2]; + + # Reshape from 1D to 2D + lo = lo.reshape(1, lo.size) + hi = hi.reshape(1, hi.size) return lo, hi Index: trunk/src/src/cshift.py =================================================================== diff --git a/trunk/src/src/cshift.py b/trunk/src/src/cshift.py --- a/trunk/src/src/cshift.py (revision 19) +++ b/trunk/src/src/cshift.py (revision 20) @@ -14,16 +14,21 @@ # y = cshift(x, m) # INPUT: # x - N-point vector - # m - amount of shift + # m - amount of shift (pos=left, neg=right) # OUTPUT: # y - vector x will be shifted by m samples to the left # # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY # http://taco.poly.edu/WaveletSoftware/ + N = x.size; - n = np.arange(N-1); + n = np.arange(N); n = np.mod(n-m, N); - y = x[n]; + + print x.shape + + y = x[0,n]; + return y \ No newline at end of file Index: trunk/src/src/deb4_basis.py =================================================================== diff --git a/trunk/src/src/deb4_basis.py b/trunk/src/src/deb4_basis.py --- a/trunk/src/src/deb4_basis.py (revision 19) +++ b/trunk/src/src/deb4_basis.py (revision 20) @@ -5,8 +5,10 @@ ''' import numpy as np -import FSfarras -import dualfilt1 +from FSfarras import * +from dualfilt1 import * +from dualtree import * +from idualtree import * def deb4_basis(N): @@ -14,25 +16,26 @@ idx = 1; J = 4; - [Faf, Fsf] = FSfarras; - [af, sf] = dualfilt1; + [Faf, Fsf] = FSfarras(); + [af, sf] = dualfilt1(); # compute transform of zero vector x = np.zeros(shape=(1,N)); - #w = dualtree(x, J, Faf, af); -# # -# # # Uses both real and imaginary wavelets -# # for i in range (1, J+1): -# # for j in range (1, 2): -# # for k in range (1, (w[i][j]).size): -# # w[i][j](k) = 1; -# # y = idualtree(w, J, Fsf, sf); -# # w[i][j](k) = 0; -# # # store it -# # Psi(:,idx) = y.T.conj(); -# # idx = idx + 1; -# # -# # # Add uniform vector (seems to be useful if there's a background -# # Psi(:,2*N+1) = 1/np.sqrt(N); -# -# return Psi \ No newline at end of file + w = dualtree(x, J, Faf, af); + + + # Uses both real and imaginary wavelets + for i in range (0, J): + for j in range (0, 1): + for k in range (0, (w[i][j]).size): + w[i][j][0,k] = 1; + y = idualtree(w, J, Fsf, sf); + w[i][j][0,k] = 0; + # store it + Psi[:,idx] = y.T.conj(); + idx = idx + 1; + + # Add uniform vector (seems to be useful if there's a background + Psi[:,2*N+1] = 1/np.sqrt(N); + + return Psi \ No newline at end of file Index: trunk/src/src/dualtree.py =================================================================== diff --git a/trunk/src/src/dualtree.py b/trunk/src/src/dualtree.py --- a/trunk/src/src/dualtree.py (revision 19) +++ b/trunk/src/src/dualtree.py (revision 20) @@ -5,7 +5,7 @@ ''' import numpy as np -import afb +from afb import * def dualtree(x, J, Faf, af): @@ -38,26 +38,39 @@ # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY # http://taco.poly.edu/WaveletSoftware/ + # ---------Trees Structure---------------# + # w [ 0 1 2 .... J ] # + # | | | | # + # [0 1] [0 1] [0 1] [0 1] # + #----------------------------------------# + # normalization x = x/np.sqrt(2); + - w = np.zeros(shape=(J,2)) ### VERIFY THIS DEFINITION + w = np.zeros(shape=(J+1), dtype=object) + for j in range (0, w.size): + w[j] = np.zeros(shape=(J+1), dtype=object) + # Tree 1 - [x1,w[1,0]] = afb(x, Faf[0,1]); # w{1}{1} + [x1, w[0][0]] = afb(x, Faf[0,0]); # w{1}{1} - for j in range (2,J): - [x1,w[j,0]] = afb(x1, af[0,1]); #check this - - w[J+1,1] = x1; + for j in range (1,J): + [x1,w[j][0]] = afb(x1, af[0,0]); ### or 0,1???? + + + + w[J][0] = x1; + # Tree 2 - [x2,w[1,2]] = afb(x, Faf[0,1]); - - for j in range (2,J): - [x2,w[j,2]] = afb(x2, af[0,1]); - - w[J+1,2] = x2; + [x2,w[0][1]] = afb(x, Faf[0,1]); + + for j in range (1,J): + [x2,w[j][1]] = afb(x2, af[0,1]); + + w[J][1] = x2; return w Index: trunk/src/src/idualtree.py =================================================================== diff --git a/trunk/src/src/idualtree.py b/trunk/src/src/idualtree.py new file mode 10644 --- /dev/null (revision 0) +++ b/trunk/src/src/idualtree.py (revision 20) @@ -0,0 +1,46 @@ +''' +Created on Jun 5, 2014 + +@author: Yolian Amaro +''' + +from sfb import * + +def idualtree(w, J, Fsf, sf): + + # Inverse Dual-tree Complex DWT + # + # USAGE: + # y = idualtree(w, J, Fsf, sf) + # INPUT: + # w - DWT coefficients + # J - number of stages + # Fsf - synthesis filters for the last stage + # sf - synthesis filters for preceeding stages + # OUTUT: + # y - output signal + # See dualtree + # + # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY + # http://taco.poly.edu/WaveletSoftware/ + + # Tree 1 + y1 = w[J][0]; + + for j in range (J-1, 0, -1): + y1 = sfb(y1, w[j][0], sf[0,0]); + + y1 = sfb(y1, w[0][0], Fsf[0,0]); + + # Tree 2 + y2 = w[J][1]; + + for j in range (J-1, 0, -1): + y2 = sfb(y2, w[j][2], sf[0,1]); + + y2 = sfb(y2, w[0][1], Fsf[0,1]); + + # normalization + y = (y1 + y2)/np.sqrt(2); + + return y Index: trunk/src/src/sfb.py =================================================================== diff --git a/trunk/src/src/sfb.py b/trunk/src/src/sfb.py new file mode 10644 --- /dev/null (revision 0) +++ b/trunk/src/src/sfb.py (revision 20) @@ -0,0 +1,68 @@ +''' +Created on Jun 5, 2014 + +@author: Yolian Amaro +''' + +from multirate import * +import numpy as np +from cshift import * + +def sfb(lo, hi, sf): + + # Synthesis filter bank + # + # USAGE: + # y = sfb(lo, hi, sf) + # INPUT: + # lo - low frqeuency input + # hi - high frequency input + # sf - synthesis filters + # sf(:, 1) - lowpass filter (even length) + # sf(:, 2) - highpass filter (even length) + # OUTPUT: + # y - output signal + # See also afb + # + # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY + # http://taco.poly.edu/WaveletSoftware/ + + N = 2*lo.size; + L = sf.size/2; + #print 'N', N + #print 'sf', sf + + + #print 'sf[:,0]', sf[:,0].shape + #print 'sf[:,1]', sf[:,1].shape + #print 'sbf hi', hi.shape + + + + # Need to change format for upfirdn funct: + lo = lo.T.conj() + lo = lo.reshape(lo.size) + + print 'sfb hi', hi + + # Need to change format for upfirdn funct: + hi = hi.T.conj() + hi = hi.reshape(hi.size) + + #hi = hi.reshape(1, hi.size) + + lo = upfirdn(lo, sf[:,0], 2, 1); + hi = upfirdn(hi, sf[:,1], 2, 1); + y = lo + hi; + y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT + y = y[0:N]; + + print 'y en sbf\n', y.shape + + y = y.reshape(1, y.size) + print 'y en sbf\n', y.shape + + y = cshift(y, 1-L/2); + + return y; +