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

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

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

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

             '''
             Created on May 26, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
   No newline at end of file
             import pywt
             import numpy as np
   No newline at end of file
             def FSfarras():
                 #function [af, sf] = FSfarras
   No newline at end of file
                 # 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
   No newline at end of file
   No newline at end of file
                 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]
                      ] );
                      ] );
  No newline at end of file
                 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]
                     ]);
                     ]);
  No newline at end of file
              No newline at end of file
-                #print a2.shape
  No newline at end of file
                 af = np.array([  [a1,a2]  ], dtype=object)
                 af = np.array([  [a1,a2]  ], dtype=object)
  No newline at end of file
                 s1 = a1[::-1]
                 s2 = a2[::-1]
                 s2 = a2[::-1]
  No newline at end of file
                 sf = np.array([  [s1,s2]  ], dtype=object)
                 sf = np.array([  [s1,s2]  ], dtype=object)
  No newline at end of file
                 return af, sf
                 return af, sf
  No newline at end of file

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

             #!/usr/bin/env python
   No newline at end of file
             #----------------------------------------------------------
             # Original MATLAB code developed by Brian Harding
             # Rewritten in python by Yolian Amaro
             # Python version 2.7
             # May 15, 2014
             # Jicamarca Radio Observatory
             #----------------------------------------------------------
   No newline at end of file
             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 *
   No newline at end of file
   No newline at end of file
             ## Calculate Forward Model
             lambda1 = 6.0
             k = 2*math.pi/lambda1
   No newline at end of file
             ## Calculate Magnetic Declination
   No newline at end of file
             #     [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012);  check this
   No newline at end of file
             # or calculate it with the above function
             dec = -1.24
   No newline at end of file
             # 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]] )
   No newline at end of file
             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]] )
   No newline at end of file
             plt.figure(1)
             plt.plot(rx, ry, 'ro')
             plt.draw()
   No newline at end of file
             # Jicamarca is nominally at a 45 degree angle
             theta = 45  - dec;
   No newline at end of file
             # 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) )
   No newline at end of file
             # Rotation matrix from antenna coord to magnetic coord (East North)
             R = np.array( [[val1, val3], [val2, val4]] );
   No newline at end of file
             # Rotate antenna positions to magnetic coord.
             AR = np.dot(R.T, antpos);
   No newline at end of file
             # Only take the East component
             r = AR[0,:]
             r.sort() # ROW VECTOR?
   No newline at end of file
             # 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?????????????????????????????
   No newline at end of file
             # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
             # background constant b.
   No newline at end of file
             # 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
   No newline at end of file
             # 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
   No newline at end of file
             # Single Gaussian
             a = np.array( [3] );
             mu = np.array( [-3.0/180*math.pi] )
             sig = np.array( [2.0/180*math.pi] )
             b = 0;
   No newline at end of file
             fact = np.zeros(shape=(Nt,1));
             factr = np.zeros(shape=(Nr,1));
   No newline at end of file
             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;
   No newline at end of file
             # # 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)));
   No newline at end of file
   No newline at end of file
             #     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)]);
   No newline at end of file
             #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
   No newline at end of file
   No newline at end of file
             ##
             # Control the type and number of inversions with:
             # SNRdBvec: the SNRs that will be used.
             # NN: the number of trials for each SNR
   No newline at end of file
             #SNRdBvec = np.linspace(5,20,10);
             SNRdBvec = np.array([15]);
             NN = 1; # number of trial at each SNR
   No newline at end of file
             # 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)
   No newline at end of file
             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);
   No newline at end of file
                     # 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);
   No newline at end of file
                     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])
   No newline at end of file
                     # 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.
   No newline at end of file
                     # This is a way of adding noise while maintaining the
                     # positive-semi-definiteness of the matrix.
   No newline at end of file
                     U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
   No newline at end of file
                     sigma_noise = (np.linalg.norm(U,'fro')/SNR);
   No newline at end of file
                     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)
   No newline at end of file
                     nz = np.multiply(temp4, temp3)
   No newline at end of file
                     #---------------------- 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))));
                     #--------------------------------------------------------------------------------
   No newline at end of file
                     Unz = U + nz;
                     Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
                     plt.figure(3);
                     plt.pcolor(abs(Rnz));
                     plt.colorbar();
   No newline at end of file
                     # Fourier Inversion ###################
                     f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
   No newline at end of file
                     for i in range(0, thetar.size):
                         th = thetar[i];
                         w = np.exp(1j*k*np.dot(r,np.sin(th)));
   No newline at end of file
                         temp = np.dot(w.T.conj(),U)
   No newline at end of file
                         f_fourier[i] = np.dot(temp, w);
   No newline at end of file
                     f_fourier = f_fourier.real; # get rid of numerical imaginary noise
   No newline at end of file
                     #print f_fourier
   No newline at end of file
   No newline at end of file
                     # Capon Inversion ######################
   No newline at end of file
                     f_capon = np.zeros(shape=(Nr,1));
   No newline at end of file
                     tic_capon = time.time();
   No newline at end of file
                     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)
   No newline at end of file
   No newline at end of file
                     toc_capon = time.time()
   No newline at end of file
                     elapsed_time_capon = toc_capon - tic_capon;
   No newline at end of file
                     f_capon = f_capon.real; # get rid of numerical imaginary noise
   No newline at end of file
                     # MaxEnt Inversion #####################
   No newline at end of file
                     # 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
   No newline at end of file
                     # 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
   No newline at end of file
                     # 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
   No newline at end of file
                     # 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;
   No newline at end of file
                     # inversion
                     F = Nr/Nt; # normalization
                     sigma = 1; # set to 1 because the difference is accounted for in G
   No newline at end of file
                     ##### 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
   No newline at end of file
                     tic_maxent = time.time();
   No newline at end of file
                     lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
   No newline at end of file
                     toc_maxent = time.time()
                     elapsed_time_maxent = toc_maxent - tic_maxent;
   No newline at end of file
                     # Whitened solution
                     def myfun(lambda1):
                         return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
   No newline at end of file
                     tic_maxEnt = time.time();
   No newline at end of file
                     #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
                     lambda1 = root(myfun,lambda0, method='krylov',  tol=1e-14);
   No newline at end of file
                     #print lambda1
                     #print lambda1.x
   No newline at end of file
                     lambda1 = lambda1.x;
   No newline at end of file
                     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);
   No newline at end of file
                     # CS inversion using Iteratively Reweighted Least Squares (IRLS)-------------
  No newline at end of file
-                    # CS inversion using irls ########################
  No newline at end of file
   No newline at end of file
                     # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
                     Psi = deb4_basis(Nr); ###### REPLACED BY LINEs BELOW (?)
-                    Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
  No newline at end of file
+             No newline at end of file
              No newline at end of file
+                    print 'FINALLY!'
              No newline at end of file
+                    print Psi.shape
  No newline at end of file
   No newline at end of file
                     # REMOVE THIS?--------------------------------
                     #wavelet1 = pywt.Wavelet('db4')
                     #Phi, Psi, x = wavelet1.wavefun(level=3)
                     # --------------------------------------------
   No newline at end of file
                     # 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 );
+             No newline at end of file
                     #         N_temp = np.array([[Nr/Nt]]);
-                    N_temp = np.array([[Nr/Nt]]);
+             No newline at end of file
                     #         g2 = np.concatenate( (gnz, N_temp), axis=0 );
-                    g2 = np.concatenate( (gnz, N_temp), axis=0 );
+             No newline at end of file
                     #         H2 = H2.T.conj();
-                    H2 = H2.T.conj();
+             No newline at end of file
              No newline at end of file
              No newline at end of file
                     #         print 'H2 shape', H2.shape
-                    print 'H2 shape', H2.shape
+             No newline at end of file
                     #         print 'Psi shape', Psi.shape
-                    print 'Psi shape', Psi.shape
+             No newline at end of file
              No newline at end of file
              No newline at end of file
                     #         s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
  No newline at end of file
-                    s = irls_dn2(H2*Psi,g2,0.5,G);
  No newline at end of file
                     #         f_cs = Psi*s;
                     #
                     #         # plot
                     #         plot(thetar,f_cs,'r.-');
                     #         hold on;
                     #         plot(thetat,fact,'k-');
                     #         hold off;
   No newline at end of file
   No newline at end of file
                     # # # Scaling and shifting
                     # # # Only necessary for capon solution
   No newline at end of file
   No newline at end of file
                     f_capon = f_capon/np.max(f_capon)*np.max(fact);
   No newline at end of file
   No newline at end of file
                     ### 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));
   No newline at end of file
   No newline at end of file
                     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);
   No newline at end of file
                     # To be able to perform dot product (align matrices) done below within the dot calculations
   No newline at end of file
   No newline at end of file
                     #f_fourier = f_fourier.T.conj()
                     #f_capon = f_capon.T.conj()
                     #f_maxent = f_maxent.T.conj()
   No newline at end of file
                     #factr = factr.T.conj()
   No newline at end of file
                     # calculate correlation
   No newline at end of file
                     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));
   No newline at end of file
   No newline at end of file
                     # calculate centered correlation
                     f0 = factr - np.mean(factr);
                     f1 = f_fourier - np.mean(f_fourier);
   No newline at end of file
                     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));
   No newline at end of file
   No newline at end of file
   No newline at end of file
                     # # # plot stuff #########################
   No newline at end of file
                     #---- 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')
   No newline at end of file
                     # 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)
   No newline at end of file
   No newline at end of file
                     #---- 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')
   No newline at end of file
                     # 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)
   No newline at end of file
                     plt.show()
   No newline at end of file
   No newline at end of file
                     # #     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);
   No newline at end of file
   No newline at end of file
                     # # 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;
   No newline at end of file
                     rmse[0,snri,Ni] = relrmse_fourier;
                     rmse[1,snri,Ni] = relrmse_capon;
                     rmse[2,snri,Ni] = relrmse_maxent;
                     #rmse[3,snri,Ni] = relrmse_cs;
   No newline at end of file
                     corrc[0,snri,Ni] = corrc_fourier;
                     corrc[1,snri,Ni] = corrc_capon;
                     corrc[2,snri,Ni] = corrc_maxent;
                     #corrc[3,snri,Ni] = corrc_cs;
   No newline at end of file
   No newline at end of file
                     print 'Capon:\t', elapsed_time_capon, 'sec';
                     print 'Maxent:\t',elapsed_time_maxent, 'sec';
                     #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
   No newline at end of file
                     print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
   No newline at end of file
                     print corr
   No newline at end of file

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

             '''
             Created on May 29, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
   No newline at end of file
             #from sp import multirate
             from cshift import *
  No newline at end of file
-            import cshift
  No newline at end of file
             from multirate import upfirdn
   No newline at end of file
             def afb(x, af):
   No newline at end of file
                 # 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/
   No newline at end of file
                 N = x.size;
                 L = (af).size/4; #L should be = 5
-                L = (af).size/2;
+             No newline at end of file
                 #print af
-                x = cshift(x,-L);
  No newline at end of file
+             No newline at end of file
+                #print 'L', L
              No newline at end of file
+                x = cshift(x,-(L-1));
              No newline at end of file
              No newline at end of file
+                #     print 'afb x', x.shape
              No newline at end of file
+                #     print 'af[:,0]',af[:,0].shape
              No newline at end of file
+                #     print 'af[:,1]',af[:,1].shape
              No newline at end of file
+                #     print '-----------------------'
  No newline at end of file
   No newline at end of file
                 # lowpass filter
                 lo = upfirdn(x, af[:,0], 1, 2);
   No newline at end of file
                 # VERIFY THIS!!!!!!!!!!!!
  No newline at end of file
                 # VERIFY THIS!!!!!!!!!!!!
                 for i in range(0, L):
  No newline at end of file
-                    lo[i] = lo[N/2+[i]] + lo[i];
  No newline at end of file
                     lo[i] = lo[N/2+i] + lo[i];
              No newline at end of file
-                lo = lo[1:N/2];
  No newline at end of file
+             No newline at end of file
+                lo = lo[0:N/2];
  No newline at end of file
   No newline at end of file
                 # highpass filter
                 # highpass filter
  No newline at end of file
-                hi = upfirdn(x, af[:,2], 1, 2);
  No newline at end of file
                 hi = upfirdn(x, af[:,1], 1, 2);
                 for j in range(0, L):
                 for j in range(0, L):
  No newline at end of file
-                    hi[j] = hi(N/2+[j]) + hi[j];
  No newline at end of file
                     hi[j] = hi[N/2+j] + hi[j];
              No newline at end of file
-                hi = hi[1:N/2];
  No newline at end of file
+             No newline at end of file
+                hi = hi[0:N/2];
              No newline at end of file
              No newline at end of file
+                # Reshape from 1D to 2D
              No newline at end of file
+                lo = lo.reshape(1, lo.size)
  No newline at end of file
                 hi = hi.reshape(1, hi.size)
                 return lo, hi
  No newline at end of file
                 return lo, hi

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

             '''
             Created on May 30, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
   No newline at end of file
             import numpy as np
   No newline at end of file
             def cshift(x, m):
   No newline at end of file
                 # Circular Shift
                 #
                 # USAGE:
                 #    y = cshift(x, m)
                 # INPUT:
                 #    x - N-point vector
                 #    m - amount of shift (pos=left, neg=right)
  No newline at end of file
-                #    m - amount of shift
  No newline at end of file
                 # 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/
   No newline at end of file
                 N = x.size;
                 N = x.size;
  No newline at end of file
-                n = np.arange(N-1);
  No newline at end of file
                 n = np.arange(N);
                 n = np.mod(n-m, N);
-                y = x[n];
  No newline at end of file
+             No newline at end of file
              No newline at end of file
+                print x.shape
              No newline at end of file
              No newline at end of file
+                y = x[0,n];
  No newline at end of file
   No newline at end of file
                 return y
  No newline at end of file
                 return y

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

             '''
             Created on May 26, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
   No newline at end of file
             import numpy as np
             from FSfarras import *
-            import FSfarras
+             No newline at end of file
             from dualfilt1 import *
-            import dualfilt1
  No newline at end of file
+             No newline at end of file
+            from dualtree import *
              No newline at end of file
+            from idualtree import *
  No newline at end of file
   No newline at end of file
             def deb4_basis(N):
   No newline at end of file
                 Psi = np.zeros(shape=(N,2*N+1));
                 idx = 1;
   No newline at end of file
                 J = 4;
                 [Faf, Fsf] = FSfarras();
-                [Faf, Fsf] = FSfarras;
+             No newline at end of file
                 [af, sf] = dualfilt1();
  No newline at end of file
-                [af, sf] = dualfilt1;
  No newline at end of file
   No newline at end of file
                 # compute transform of zero vector
                 x = np.zeros(shape=(1,N));
                 w = dualtree(x, J, Faf, af);
-                #w = dualtree(x, J, Faf, af);
+             No newline at end of file
              No newline at end of file
-            #     #
+             No newline at end of file
              No newline at end of file
-            #     #     # Uses both real and imaginary wavelets
+             No newline at end of file
                 # Uses both real and imaginary wavelets
-            #     #     for i in range (1, J+1):
+             No newline at end of file
                 for i in range (0, J):
-            #     #         for j in range (1, 2):
+             No newline at end of file
                     for j in range (0, 1):
-            #     #             for k in range (1, (w[i][j]).size):
+             No newline at end of file
                         for k in range (0, (w[i][j]).size):
-            #     #                 w[i][j](k) = 1;
+             No newline at end of file
                             w[i][j][0,k] = 1;
-            #     #                 y = idualtree(w, J, Fsf, sf);
+             No newline at end of file
                             y = idualtree(w, J, Fsf, sf);
-            #     #                 w[i][j](k) = 0;
+             No newline at end of file
                             w[i][j][0,k] = 0;
-            #     #                 # store it
+             No newline at end of file
                             # store it
-            #     #                 Psi(:,idx) = y.T.conj();
+             No newline at end of file
                             Psi[:,idx] = y.T.conj();
-            #     #                 idx = idx + 1;
+             No newline at end of file
                             idx = idx + 1;
-            #     #
+             No newline at end of file
              No newline at end of file
-            #     #     # Add uniform vector (seems to be useful if there's a background
+             No newline at end of file
                 # Add uniform vector (seems to be useful if there's a background
-            #     #     Psi(:,2*N+1) = 1/np.sqrt(N);
+             No newline at end of file
                 Psi[:,2*N+1] = 1/np.sqrt(N);
              No newline at end of file
              No newline at end of file
-            #     return Psi
  No newline at end of file
+             No newline at end of file
+                return Psi
  No newline at end of file

trunk/src/src/dualtree.py +23 -10

             '''
             Created on May 29, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
   No newline at end of file
             import numpy as np
             from afb import *
  No newline at end of file
-            import afb
  No newline at end of file
   No newline at end of file
             def dualtree(x, J, Faf, af):
   No newline at end of file
                 # 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/
   No newline at end of file
+                # ---------Trees Structure---------------#
              No newline at end of file
+                # w [ 0      1      2     ....    J ]    #
              No newline at end of file
+                #     |      |      |             |      #
              No newline at end of file
+                #   [0 1]  [0 1]  [0 1]         [0 1]    #
              No newline at end of file
+                #----------------------------------------#
              No newline at end of file
+  No newline at end of file
                 # normalization
                 x = x/np.sqrt(2);
              No newline at end of file
-                w = np.zeros(shape=(J,2)) ### VERIFY THIS DEFINITION
  No newline at end of file
+             No newline at end of file
+                w = np.zeros(shape=(J+1), dtype=object)
              No newline at end of file
              No newline at end of file
+                for j in range (0, w.size):
              No newline at end of file
+                    w[j] = np.zeros(shape=(J+1), dtype=object)
  No newline at end of file
   No newline at end of file
                 # Tree 1
                 [x1, w[0][0]] = afb(x, Faf[0,0]); # w{1}{1}
  No newline at end of file
-                [x1,w[1,0]] = afb(x, Faf[0,1]); # w{1}{1}
  No newline at end of file
              No newline at end of file
-                for j in range (2,J):
              No newline at end of file
-                    [x1,w[j,0]] = afb(x1, af[0,1]); #check this
  No newline at end of file
                 for j in range (1,J):
-                w[J+1,1] = x1;
  No newline at end of file
+             No newline at end of file
+                    [x1,w[j][0]] = afb(x1, af[0,0]);  ### or 0,1????
              No newline at end of file
              No newline at end of file
              No newline at end of file
              No newline at end of file
+                w[J][0] = x1;
  No newline at end of file
   No newline at end of file
                 # Tree 2
                 [x2,w[0][1]] = afb(x, Faf[0,1]);
  No newline at end of file
-                [x2,w[1,2]] = afb(x, Faf[0,1]);
  No newline at end of file
                 for j in range (1,J):
-                for j in range (2,J):
+             No newline at end of file
                     [x2,w[j][1]] = afb(x2, af[0,1]);
  No newline at end of file
-                    [x2,w[j,2]] = afb(x2, af[0,1]);
  No newline at end of file
                 w[J][1] = x2;
  No newline at end of file
-                w[J+1,2] = x2;
  No newline at end of file
   No newline at end of file
                 return w
   No newline at end of file

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