jro_soft/compressed_sensing Commit - r19:20

Fixed errors, implemented sfb and idualtree classes

yamaro -

r19:20

parent child

trunk/src/src/idualtree.py created 10644 +46 0

@@ -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

trunk/src/src/sfb.py created 10644 +68 0

@@ -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

trunk/src/src/FSfarras.py 0 -2

             '''
             Created on May 26, 2014
             @author: Yolian Amaro
             '''
             import pywt
             import numpy as np
             def FSfarras():
                 #function [af, sf] = FSfarras
                 # Farras filters organized for the dual-tree
                 # complex DWT.
                 #
                 # USAGE:
                 #    [af, sf] = FSfarras
                 # OUTPUT:
                 #    af{i}, i = 1,2 - analysis filters for tree i
                 #    sf{i}, i = 1,2 - synthesis filters for tree i
                 # See farras, dualtree, dualfilt1.
                 #
                 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
                 # http://taco.poly.edu/WaveletSoftware/
                 #
                 # Translated to Python by Yolian Amaro
                 a1 = np.array( [
                       [                0,                  0],
                       [-0.08838834764832,   -0.01122679215254],
                       [ 0.08838834764832,   0.01122679215254],
                       [ 0.69587998903400,    0.08838834764832],
                       [ 0.69587998903400,   0.08838834764832],
                       [ 0.08838834764832,   -0.69587998903400],
                       [-0.08838834764832,   0.69587998903400],
                       [ 0.01122679215254,   -0.08838834764832],
                       [ 0.01122679215254,   -0.08838834764832],
                       [0,                                   0]
                      ] );
                 a2 = np.array([
                       [ 0.01122679215254,                  0],
                       [ 0.01122679215254,                  0],
                       [-0.08838834764832,  -0.08838834764832],
                       [ 0.08838834764832,  -0.08838834764832],
                       [ 0.69587998903400,   0.69587998903400],
                       [ 0.69587998903400,  -0.69587998903400],
                       [ 0.08838834764832,   0.08838834764832],
                       [-0.08838834764832,   0.08838834764832],
                       [                0,   0.01122679215254],
                       [                0,  -0.01122679215254]
                     ]);
-                #print a2.shape
                 af = np.array([  [a1,a2]  ], dtype=object)
                 s1 = a1[::-1]
                 s2 = a2[::-1]
                 sf = np.array([  [s1,s2]  ], dtype=object)
                 return af, sf

trunk/src/src/RICS_paper_compare_noise.py +14 -11

             #!/usr/bin/env python
             #----------------------------------------------------------
             # Original MATLAB code developed by Brian Harding
             # Rewritten in python by Yolian Amaro
             # Python version 2.7
             # May 15, 2014
             # Jicamarca Radio Observatory
             #----------------------------------------------------------
             import math
             import numpy as np
             import matplotlib.pyplot as plt
             from scipy import linalg
             import time
             from y_hysell96 import*
             from deb4_basis import *
             from modelf import *
             #from scipy.optimize import fsolve
             from scipy.optimize import root
             import pywt
             from irls_dn2 import *
             ## Calculate Forward Model
             lambda1 = 6.0
             k = 2*math.pi/lambda1
             ## Calculate Magnetic Declination
             #     [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012);  check this
             # or calculate it with the above function
             dec = -1.24
             # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
             rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
             ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
             antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240],
                             [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
             plt.figure(1)
             plt.plot(rx, ry, 'ro')
             plt.draw()
             # Jicamarca is nominally at a 45 degree angle
             theta = 45  - dec;
             # Rotation matrix from antenna coord to magnetic coord (East North)
             theta_rad = math.radians(theta)        #  trig functions take radians as argument
             val1 = float( math.cos(theta_rad) )
             val2 = float( math.sin(theta_rad) )
             val3 = float( -1*math.sin(theta_rad))
             val4 = float( math.cos(theta_rad) )
             # Rotation matrix from antenna coord to magnetic coord (East North)
             R = np.array( [[val1, val3], [val2, val4]] );
             # Rotate antenna positions to magnetic coord.
             AR = np.dot(R.T, antpos);
             # Only take the East component
             r = AR[0,:]
             r.sort() # ROW VECTOR?
             # Truth model (high and low resolution)
             Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
             thbound = 9.0/180*math.pi; # the width of the domain in angle space
             thetat = np.linspace(-thbound, thbound,Nt) # image domain
             thetat = np.transpose(thetat) # transpose  # FUNCIONA??????????????????????????????
             Nr = (256.0); # number of pixels in reconstructed image: low res
             thetar = np.linspace(-thbound, thbound,Nr)  # reconstruction domain
             thetar = np.transpose(thetar) #transpose   # FUNCIONA?????????????????????????????
             # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
             # background constant b.
             # Triple Gaussian
             # a = np.array([3, 5, 2]);
             # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
             # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
             # b = 0; # background
             # Double Gaussian
             # a = np.array([3, 5]);
             # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
             # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
             # b = 0; # background
             # Single Gaussian
             a = np.array( [3] );
             mu = np.array( [-3.0/180*math.pi] )
             sig = np.array( [2.0/180*math.pi] )
             b = 0;
             fact = np.zeros(shape=(Nt,1));
             factr = np.zeros(shape=(Nr,1));
             for i in range(0, a.size):
                 temp = (-(thetat-mu[i])**2/(sig[i]**2))
                 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
                 for j in range(0, temp.size):
                     fact[j] = fact[j] + a[i]*math.exp(temp[j]);
                 for m in range(0, tempr.size):
                     factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
             fact = fact + b;
             factr = factr + b;
             # # model for f: Square pulse
             # for j in range(0, fact.size):
             #     if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi):
             #         fact[j] = 0
             #     else:
             #         fact[j] = 1
             # for k in range(0, factr.size):
             #     if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
             #         fact[k] = 0
             #     else:
             #         fact[k] = 1
             #
             #
             # # model for f: triangle pulse
             # mu = -1.0/180*math.pi;
             # sig = 5.0/180*math.pi;
             # wind1 = theta > mu-sig and theta < mu;
             # wind2 = theta < mu+sig and theta > mu;
             # fact = wind1 * (theta - (mu - sig));
             # factr = wind1 * (thetar - (mu - sig));
             # fact = fact + wind2 * (-(theta-(mu+sig)));
             # factr = factr + wind2 * (-(thetar-(mu+sig)));
             #     fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1
             I = sum(fact)[0];
             fact = fact/I; # normalize to sum(f)==1
             factr = factr/I; # normalize to sum(f)==1
             #plt.figure()
             #plt.plot(thetat,fact,'r');
             #plt.plot(thetar,factr,'k.');
             #xlim([min(thetat) max(thetat)]);
             #x = np.linspace(thetat.min(), thetat.max) ????
             #for i in range(0, thetat.size):
             plt.figure(2)
             plt.plot(thetat, fact, 'r--')
             plt.plot(thetar, factr, 'ro')
             plt.draw()
             # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
             ##
             # Control the type and number of inversions with:
             # SNRdBvec: the SNRs that will be used.
             # NN: the number of trials for each SNR
             #SNRdBvec = np.linspace(5,20,10);
             SNRdBvec = np.array([15]);
             NN = 1; # number of trial at each SNR
             # if using vector arguments should be: (4,SNRdBvec.size,NN)
             corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
             corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
             rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
             for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
                 for Ni in range(0, NN):
                     SNRdB = SNRdBvec[snri];
                     SNR = 10**(SNRdB/10.0);
                     # Calculate cross-correlation matrix (Fourier components of image)
                     # This is an inefficient way to do this.
                     R = np.zeros(shape=(r.size, r.size), dtype=object);
                     for i1 in range(0, r.size):
                         for i2 in range(0,r.size):
                             R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat))))
                             R[i1,i2] = sum(R[i1,i2])
                     # Add uncertainty
                     # This is an ad-hoc way of adding "noise".  It models some combination of
                     # receiver noise and finite integration times.  We could use a more
                     # advanced model (like in Yu et al 2000) in the future.
                     # This is a way of adding noise while maintaining the
                     # positive-semi-definiteness of the matrix.
                     U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
                     sigma_noise = (np.linalg.norm(U,'fro')/SNR);
                     temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
                     temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
                     temp3 = ((abs(U) > 0).astype(float))   # upper triangle of 1's
                     temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
                     nz = np.multiply(temp4, temp3)
                     #---------------------- Eliminar esto:------------------------------------------
                     #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int))
                     #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
                     #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
                     #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))));
                     #--------------------------------------------------------------------------------
                     Unz = U + nz;
                     Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
                     plt.figure(3);
                     plt.pcolor(abs(Rnz));
                     plt.colorbar();
                     # Fourier Inversion ###################
                     f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
                     for i in range(0, thetar.size):
                         th = thetar[i];
                         w = np.exp(1j*k*np.dot(r,np.sin(th)));
                         temp = np.dot(w.T.conj(),U)
                         f_fourier[i] = np.dot(temp, w);
                     f_fourier = f_fourier.real; # get rid of numerical imaginary noise
                     #print f_fourier
                     # Capon Inversion ######################
                     f_capon = np.zeros(shape=(Nr,1));
                     tic_capon = time.time();
                     for i in range(0, thetar.size):
                         th = thetar[i];
                         w = np.exp(1j*k*np.dot(r,np.sin(th)));
                         f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
                     toc_capon = time.time()
                     elapsed_time_capon = toc_capon - tic_capon;
                     f_capon = f_capon.real; # get rid of numerical imaginary noise
                     # MaxEnt Inversion #####################
                     # create the appropriate sensing matrix (split into real and imaginary # parts)
                     M = (r.size-1)*(r.size);
                     Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
                     Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
                     # need to re-index our measurements from matrix R into vector g
                     g = np.zeros(shape=(M,1));
                     gnz = np.zeros(shape=(M,1)); # noisy version of g
                     # triangular indexing to perform this re-indexing
                     T = np.ones(shape=(r.size,r.size));
                     [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
                     # build H
                     for i1 in range(0, r.size):
                         for i2 in range(i1+1, r.size):
                             idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward
                             idx1 = 2*idx; # because index starts at 0
                             idx2 = 2*idx+1;
                             Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T;
                             Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T;
                             Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
                             Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
                             g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later
                             g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again
                             gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
                             gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
                     # inversion
                     F = Nr/Nt; # normalization
                     sigma = 1; # set to 1 because the difference is accounted for in G
                     ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
                     G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2
                     tic_maxent = time.time();
                     lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
                     toc_maxent = time.time()
                     elapsed_time_maxent = toc_maxent - tic_maxent;
                     # Whitened solution
                     def myfun(lambda1):
                         return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
                     tic_maxEnt = time.time();
                     #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
                     lambda1 = root(myfun,lambda0, method='krylov',  tol=1e-14);
                     #print lambda1
                     #print lambda1.x
                     lambda1 = lambda1.x;
                     toc_maxEnt = time.time();
                     f_maxent = modelf(lambda1, Hr, F);
                     ystar = myfun(lambda1);
                     Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G));
                     ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
                     es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
                     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')
                     #Phi, Psi, x = wavelet1.wavefun(level=3)
                     # --------------------------------------------
                     # add "sum to 1" constraint
-                    H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
+                    #         H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
-                    N_temp = np.array([[Nr/Nt]]);
+                    #         N_temp = np.array([[Nr/Nt]]);
-                    g2 = np.concatenate( (gnz, N_temp), axis=0 );
+                    #         g2 = np.concatenate( (gnz, N_temp), axis=0 );
-                    H2 = H2.T.conj();
+                    #         H2 = H2.T.conj();
+                    #
-                    print 'H2 shape', H2.shape
+                    #         print 'H2 shape', H2.shape
-                    print 'Psi shape', Psi.shape
+                    #         print 'Psi shape', Psi.shape
+                    #
-                    s = irls_dn2(H2*Psi,g2,0.5,G);
+                    #         s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
                     #         f_cs = Psi*s;
                     #
                     #         # plot
                     #         plot(thetar,f_cs,'r.-');
                     #         hold on;
                     #         plot(thetat,fact,'k-');
                     #         hold off;
                     # # # Scaling and shifting
                     # # # Only necessary for capon solution
                     f_capon = f_capon/np.max(f_capon)*np.max(fact);
                     ### analyze stuff ######################
                     # calculate MSE
                     rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
                     rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
                     rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
                     #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
                     relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
                     relrmse_capon = rmse_capon / np.linalg.norm(fact);
                     relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
                     #relrmse_cs = rmse_cs / np.norm(fact);
                     # To be able to perform dot product (align matrices) done below within the dot calculations
                     #f_fourier = f_fourier.T.conj()
                     #f_capon = f_capon.T.conj()
                     #f_maxent = f_maxent.T.conj()
                     #factr = factr.T.conj()
                     # calculate correlation
                     corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
                     corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
                     corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
                     #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
                     # calculate centered correlation
                     f0 = factr - np.mean(factr);
                     f1 = f_fourier - np.mean(f_fourier);
                     corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
                     f1 = f_capon - np.mean(f_capon);
                     corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
                     f1 = f_maxent - np.mean(f_maxent);
                     corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
                     #f1 = f_cs - mean(f_cs);
                     #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
                     # # # plot stuff #########################
                     #---- Capon----
                     plt.figure(4)
                     plt.subplot(2, 1, 1)
                     plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
                     plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
                     plt.ylabel('Power (arbitrary units)')
                     plt.legend(loc='upper right')
                     # formatting y-axis
                     locs,labels = plt.yticks()
                     plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
                     #---- MaxEnt----
                     plt.subplot(2, 1, 2)
                     plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
                     plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
                     plt.ylabel('Power (arbitrary units)')
                     plt.legend(loc='upper right')
                     # formatting y-axis
                     locs,labels = plt.yticks()
                     plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
                     plt.show()
                     # #     PLOT PARA COMPRESSED SENSING
                     # #
                     # #     subplot(3,1,3);
                     # #         plot(180/pi*thetar,f_cs,'r-');
                     # #         hold on;
                     # #         plot(180/pi*thetat,fact,'k--');
                     # #         hold off;
                     # #         ylim([min(f_cs) 1.1*max(fact)]);
                     # # #         title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
                     # # #         title 'Compressed Sensing - Debauchies Wavelets'
                     # #         xlabel 'Degrees'
                     # #         ylabel({'Power';'(arbitrary units)'})
                     # #         legend('Comp. Sens.','Truth');
                     # #
                     # # #     set(gcf,'Position',[749   143   528   881]); # CSL
                     # # #     set(gcf,'Position',[885   -21   528   673]); # macbook
                     # #     pause(0.01);
                     # # Store Results
                     corr[0, snri, Ni] = corr_fourier;
                     corr[1, snri, Ni] = corr_capon;
                     corr[2, snri, Ni] = corr_maxent;
                     #corr[3, snri, Ni] = corr_cs;
                     rmse[0,snri,Ni] = relrmse_fourier;
                     rmse[1,snri,Ni] = relrmse_capon;
                     rmse[2,snri,Ni] = relrmse_maxent;
                     #rmse[3,snri,Ni] = relrmse_cs;
                     corrc[0,snri,Ni] = corrc_fourier;
                     corrc[1,snri,Ni] = corrc_capon;
                     corrc[2,snri,Ni] = corrc_maxent;
                     #corrc[3,snri,Ni] = corrc_cs;
                     print 'Capon:\t', elapsed_time_capon, 'sec';
                     print 'Maxent:\t',elapsed_time_maxent, 'sec';
                     #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
                     print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
                     print corr
   No newline at end of file

trunk/src/src/afb.py +21 -8

             '''
             Created on May 29, 2014
             @author: Yolian Amaro
             '''
             #from sp import multirate
-            import cshift
+            from cshift import *
             from multirate import upfirdn
             def afb(x, af):
                 # Analysis filter bank
                 #
                 # USAGE:
                 #    [lo, hi] = afb(x, af)
                 # INPUT:
                 #    x - N-point vector, where
                 #            1) N is even
                 #            2) N >= length(af)
                 #    af - analysis filters
                 #    af(:, 1) - lowpass filter (even length)
                 #    af(:, 2) - highpass filter (even length)
                 # OUTPUT:
                 #    lo - Low frequecy output
                 #    hi - High frequency output
                 # EXAMPLE:
                 #    [af, sf] = farras;
                 #    x = rand(1,64);
                 #    [lo, hi] = afb(x, af);
                 #    y = sfb(lo, hi, sf);
                 #    err = x - y;
                 #    max(abs(err))
                 #
                 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
                 # http://taco.poly.edu/WaveletSoftware/
                 N = x.size;
-                L = (af).size/2;
+                L = (af).size/4; #L should be = 5
-                x = cshift(x,-L);
+                #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

trunk/src/src/cshift.py +8 -3

             '''
             Created on May 30, 2014
             @author: Yolian Amaro
             '''
             import numpy as np
             def cshift(x, m):
                 # Circular Shift
                 #
                 # USAGE:
                 #    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

trunk/src/src/deb4_basis.py +24 -21

             '''
             Created on May 26, 2014
             @author: Yolian Amaro
             '''
             import numpy as np
-            import FSfarras
+            from FSfarras import *
-            import dualfilt1
+            from dualfilt1 import *
+            from dualtree import *
+            from idualtree import *
             def deb4_basis(N):
                 Psi = np.zeros(shape=(N,2*N+1));
                 idx = 1;
                 J = 4;
-                [Faf, Fsf] = FSfarras;
+                [Faf, Fsf] = FSfarras();
-                [af, sf] = dualfilt1;
+                [af, sf] = dualfilt1();
                 # compute transform of zero vector
                 x = np.zeros(shape=(1,N));
-                #w = dualtree(x, J, Faf, af);
+                w = dualtree(x, J, Faf, af);
-            #     #
-            #     #     # Uses both real and imaginary wavelets
-            #     #     for i in range (1, J+1):
+                # Uses both real and imaginary wavelets
-            #     #         for j in range (1, 2):
+                for i in range (0, J):
-            #     #             for k in range (1, (w[i][j]).size):
+                    for j in range (0, 1):
-            #     #                 w[i][j](k) = 1;
+                        for k in range (0, (w[i][j]).size):
-            #     #                 y = idualtree(w, J, Fsf, sf);
+                            w[i][j][0,k] = 1;
-            #     #                 w[i][j](k) = 0;
+                            y = idualtree(w, J, Fsf, sf);
-            #     #                 # store it
+                            w[i][j][0,k] = 0;
-            #     #                 Psi(:,idx) = y.T.conj();
+                            # store it
-            #     #                 idx = idx + 1;
+                            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);
+                # 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
+                return Psi
  No newline at end of file

trunk/src/src/dualtree.py +26 -13

             '''
             Created on May 29, 2014
             @author: Yolian Amaro
             '''
             import numpy as np
-            import afb
+            from afb import *
             def dualtree(x, J, Faf, af):
                 # Dual-tree Complex Discrete Wavelet Transform
                 #
                 # USAGE:
                 #    w = dualtree(x, J, Faf, af)
                 # INPUT:
                 #   x - N-point vector
                 #       1) N is divisible by 2^J
                 #       2) N >= 2^(J-1)*length(af)
                 #   J - number of stages
                 #   Faf - filters for the first stage
                 #   af - filters for the remaining stages
                 # OUTPUT:
                 #   w - DWT coefficients
                 #      w{j}{1}, j = 1..J - real part
                 #      w{j}{2}, j = 1..J - imaginary part
                 #      w{J+1}{d} - lowpass coefficients, d = 1,2
                 # EXAMPLE:
                 #    x = rand(1, 512);
                 #    J = 4;
                 #    [Faf, Fsf] = FSfarras;
                 #    [af, sf] = dualfilt1;
                 #    w = dualtree(x, J, Faf, af);
                 #    y = idualtree(w, J, Fsf, sf);
                 #    err = x - y;
                 #    max(abs(err))
                 #
                 # 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]);
+                [x2,w[0][1]] = afb(x, Faf[0,1]);
-                for j in range (2,J):
+                for j in range (1,J):
-                    [x2,w[j,2]] = afb(x2, af[0,1]);
+                    [x2,w[j][1]] = afb(x2, af[0,1]);
-                w[J+1,2] = x2;
+                w[J][1] = x2;
                 return w

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages