nikcleju@9: # -*- coding: utf-8 -*- nikcleju@9: """ nikcleju@9: Created on Sat Nov 05 18:39:54 2011 nikcleju@9: nikcleju@9: @author: Nic nikcleju@9: """ nikcleju@9: nikcleju@9: import numpy as np nikcleju@9: nikcleju@9: #function s=SL0(A, x, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s) nikcleju@9: def SL0(A, x, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, A_pinv=None, true_s=None): nikcleju@9: # nikcleju@9: # SL0(A, x, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s) nikcleju@9: # nikcleju@9: # Returns the sparsest vector s which satisfies underdetermined system of nikcleju@9: # linear equations A*s=x, using Smoothed L0 (SL0) algorithm. Note that nikcleju@9: # the matrix A should be a 'wide' matrix (more columns than rows). The nikcleju@9: # number of the rows of matrix A should be equal to the length of the nikcleju@9: # column vector x. nikcleju@9: # nikcleju@9: # The first 3 arguments should necessarily be provided by the user. The nikcleju@9: # other parameters have defult values calculated within the function, or nikcleju@9: # may be provided by the user. nikcleju@9: # nikcleju@9: # Sequence of Sigma (sigma_min and sigma_decrease_factor): nikcleju@9: # This is a decreasing geometric sequence of positive numbers: nikcleju@9: # - The first element of the sequence of sigma is calculated nikcleju@9: # automatically. The last element is given by 'sigma_min', and the nikcleju@9: # change factor for decreasing sigma is given by 'sigma_decrease_factor'. nikcleju@9: # - The default value of 'sigma_decrease_factor' is 0.5. Larger value nikcleju@9: # gives better results for less sparse sources, but it uses more steps nikcleju@9: # on sigma to reach sigma_min, and hence it requires higher nikcleju@9: # computational cost. nikcleju@9: # - There is no default value for 'sigma_min', and it should be nikcleju@9: # provided by the user (depending on his/her estimated source noise nikcleju@9: # level, or his/her desired accuracy). By `noise' we mean here the nikcleju@9: # noise in the sources, that is, the energy of the inactive elements of nikcleju@9: # 's'. For example, by the noiseless case, we mean the inactive nikcleju@9: # elements of 's' are exactly equal to zero. As a rule of tumb, for the nikcleju@9: # noisy case, sigma_min should be about 2 to 4 times of the standard nikcleju@9: # deviation of this noise. For the noiseless case, smaller 'sigma_min' nikcleju@9: # results in better estimation of the sparsest solution, and hence its nikcleju@9: # value is determined by the desired accuracy. nikcleju@9: # nikcleju@9: # mu_0: nikcleju@9: # The value of mu_0 scales the sequence of mu. For each vlue of nikcleju@9: # sigma, the value of mu is chosen via mu=mu_0*sigma^2. Note that this nikcleju@9: # value effects Convergence. nikcleju@9: # The default value is mu_0=2 (see the paper). nikcleju@9: # nikcleju@9: # L: nikcleju@9: # number of iterations of the internal (steepest ascent) loop. The nikcleju@9: # default value is L=3. nikcleju@9: # nikcleju@9: # A_pinv: nikcleju@9: # is the pseudo-inverse of matrix A defined by A_pinv=A'*inv(A*A'). nikcleju@9: # If it is not provided, it will be calculated within the function. If nikcleju@9: # you use this function for solving x(t)=A s(t) for different values of nikcleju@9: # 't', it would be a good idea to calculate A_pinv outside the function nikcleju@9: # to prevent its re-calculation for each 't'. nikcleju@9: # nikcleju@9: # true_s: nikcleju@9: # is the true value of the sparse solution. This argument is for nikcleju@9: # simulation purposes. If it is provided by the user, then the function nikcleju@9: # will calculate the SNR of the estimation for each value of sigma and nikcleju@9: # it provides a progress report. nikcleju@9: # nikcleju@9: # Authors: Massoud Babaie-Zadeh and Hossein Mohimani nikcleju@9: # Version: 1.4 nikcleju@9: # Last modified: 4 April 2010. nikcleju@9: # nikcleju@9: # nikcleju@9: # Web-page: nikcleju@9: # ------------------ nikcleju@9: # http://ee.sharif.ir/~SLzero nikcleju@9: # nikcleju@9: # Code History: nikcleju@9: #-------------- nikcleju@9: # Version 2.0: 4 April 2010 nikcleju@9: # Doing a few small modifications that enable the code to work also nikcleju@9: # for complex numbers (not only for real numbers). nikcleju@9: # nikcleju@9: # Version 1.3: 13 Sep 2008 nikcleju@9: # Just a few modifications in the comments nikcleju@9: # nikcleju@9: # Version 1.2: Adding some more comments in the help section nikcleju@9: # nikcleju@9: # Version 1.1: 4 August 2008 nikcleju@9: # - Using MATLAB's pseudo inverse function to generalize for the case nikcleju@9: # the matrix A is not full-rank. nikcleju@9: # nikcleju@9: # Version 1.0 (first official version): 4 July 2008. nikcleju@9: # nikcleju@9: # First non-official version and algorithm development: Summer 2006 nikcleju@9: #--------------------- nikcleju@9: # Original Matab code: nikcleju@9: # nikcleju@9: #if nargin < 4 nikcleju@9: # sigma_decrease_factor = 0.5; nikcleju@9: # A_pinv = pinv(A); nikcleju@9: # mu_0 = 2; nikcleju@9: # L = 3; nikcleju@9: # ShowProgress = logical(0); nikcleju@9: #elseif nargin == 4 nikcleju@9: # A_pinv = pinv(A); nikcleju@9: # mu_0 = 2; nikcleju@9: # L = 3; nikcleju@9: # ShowProgress = logical(0); nikcleju@9: #elseif nargin == 5 nikcleju@9: # A_pinv = pinv(A); nikcleju@9: # L = 3; nikcleju@9: # ShowProgress = logical(0); nikcleju@9: #elseif nargin == 6 nikcleju@9: # A_pinv = pinv(A); nikcleju@9: # ShowProgress = logical(0); nikcleju@9: #elseif nargin == 7 nikcleju@9: # ShowProgress = logical(0); nikcleju@9: #elseif nargin == 8 nikcleju@9: # ShowProgress = logical(1); nikcleju@9: #else nikcleju@9: # error('Error in calling SL0 function'); nikcleju@9: #end nikcleju@9: # nikcleju@9: # nikcleju@9: ## Initialization nikcleju@9: ##s = A\x; nikcleju@9: #s = A_pinv*x; nikcleju@9: #sigma = 2*max(abs(s)); nikcleju@9: # nikcleju@9: ## Main Loop nikcleju@9: #while sigma>sigma_min nikcleju@9: # for i=1:L nikcleju@9: # delta = OurDelta(s,sigma); nikcleju@9: # s = s - mu_0*delta; nikcleju@9: # s = s - A_pinv*(A*s-x); # Projection nikcleju@9: # end nikcleju@9: # nikcleju@9: # if ShowProgress nikcleju@9: # fprintf(' sigma=#f, SNR=#f\n',sigma,estimate_SNR(s,true_s)) nikcleju@9: # end nikcleju@9: # nikcleju@9: # sigma = sigma * sigma_decrease_factor; nikcleju@9: #end nikcleju@9: # nikcleju@9: # nikcleju@9: ##################################################################### nikcleju@9: #function delta=OurDelta(s,sigma) nikcleju@9: # nikcleju@9: #delta = s.*exp(-abs(s).^2/sigma^2); nikcleju@9: # nikcleju@9: ##################################################################### nikcleju@9: #function SNR=estimate_SNR(estim_s,true_s) nikcleju@9: # nikcleju@9: #err = true_s - estim_s; nikcleju@9: #SNR = 10*log10(sum(abs(true_s).^2)/sum(abs(err).^2)); nikcleju@9: nikcleju@9: # End of original Matab code nikcleju@9: #---------------------------- nikcleju@9: nikcleju@9: if A_pinv is None: nikcleju@9: A_pinv = np.linalg.pinv(A) nikcleju@9: nikcleju@9: if true_s is not None: nikcleju@9: ShowProgress = True nikcleju@9: else: nikcleju@9: ShowProgress = False nikcleju@9: nikcleju@9: # Initialization nikcleju@9: #s = A\x; nikcleju@9: s = np.dot(A_pinv,x) nikcleju@9: sigma = 2.0*np.abs(s).max() nikcleju@9: nikcleju@9: # Main Loop nikcleju@9: while sigma>sigma_min: nikcleju@9: for i in np.arange(L): nikcleju@9: delta = OurDelta(s,sigma) nikcleju@9: s = s - mu_0*delta nikcleju@9: s = s - np.dot(A_pinv,(np.dot(A,s)-x)) # Projection nikcleju@9: nikcleju@9: if ShowProgress: nikcleju@9: #fprintf(' sigma=#f, SNR=#f\n',sigma,estimate_SNR(s,true_s)) nikcleju@9: string = ' sigma=%f, SNR=%f\n' % sigma,estimate_SNR(s,true_s) nikcleju@9: print string nikcleju@9: nikcleju@9: sigma = sigma * sigma_decrease_factor nikcleju@9: nikcleju@9: return s nikcleju@9: nikcleju@9: nikcleju@9: #################################################################### nikcleju@9: #function delta=OurDelta(s,sigma) nikcleju@9: def OurDelta(s,sigma): nikcleju@9: nikcleju@9: return s * np.exp( (-np.abs(s)**2) / sigma**2) nikcleju@9: nikcleju@9: #################################################################### nikcleju@9: #function SNR=estimate_SNR(estim_s,true_s) nikcleju@9: def estimate_SNR(estim_s, true_s): nikcleju@9: nikcleju@9: err = true_s - estim_s nikcleju@9: return 10*np.log10((true_s**2).sum()/(err**2).sum()) nikcleju@9: