demo_streamingDWT

 demo_streamingDWT

 Solves the following streaming problem
 min_a  \|a\|_1 + 1/2*||Phi Psi a + y||_2^2

 with Psi as an overlapping representation matrix and 
 by adding and removing one set of measurements at every time instant

 Applications:
           streaming signal recovery using (overlapping) DWT
           compare recovery using block-based DWT (sym=0) and
           overlapping DWT (sym = 3) near line 80
 We can also add any other regularization operator in the reconstruction

 Written by: Salman Asif, Georgia Tech
 Email: sasif@gatech.edu
 Created: March 2013

 TODO: this code needs some debugging (not stable)...

0001 % demo_streamingDWT
0002 %
0003 % Solves the following streaming problem
0004 % min_a  \|a\|_1 + 1/2*||Phi Psi a + y||_2^2
0005 %
0006 % with Psi as an overlapping representation matrix and
0007 % by adding and removing one set of measurements at every time instant
0008 %
0009 % Applications:
0010 %           streaming signal recovery using (overlapping) DWT
0011 %           compare recovery using block-based DWT (sym=0) and
0012 %           overlapping DWT (sym = 3) near line 80
0013 % We can also add any other regularization operator in the reconstruction
0014 %
0015 % Written by: Salman Asif, Georgia Tech
0016 % Email: sasif@gatech.edu
0017 % Created: March 2013
0018 %
0019 % TODO: this code needs some debugging (not stable)...
0020 
0021 clear
0022 close all force
0023  
0024 % Limit the number of computational threads (for profiling)
0025 % maxNumCompThreads(1);
0026 
0027 %% Setup path
0028 mname = mfilename;
0029 mpath = mfilename('fullpath');
0030 mdir = mpath(1:end-length(mname));
0031 cd(mdir);
0032 
0033 addpath utils/
0034 addpath utils/utils_Wavelet
0035 addpath utils/utils_LOT
0036 addpath solvers/
0037 
0038 fprintf(['----------',datestr(now),'-------%s------------------\n'],mname)
0039 
0040 % load RandomStates
0041 %
0042 rseed = 2013;
0043 rand('state',rseed);
0044 randn('state',rseed);
0045 
0046 % simulation parameters
0047 mType = 'sign'; % {'randn','orth','rdct','streaming','subsample'};
0048 sType = 'pcwreg'; % % menu('Choose function:', 'heavisine', 'pcwpoly (N>256,daub4)', 'pcwreg (N>128,daub4)', 'doppler', ...);
0049 SNR = 35;       % additive Gaussian noise
0050 
0051 % system parameters
0052 N = 256;   % signal length
0053 R = 4;      % compression rate
0054 M = round(N/R);    % no. of measurements
0055 
0056 LM = N; % LM: length of measurement window
0057 
0058 % streaming window
0059 P = 5; % size of the working window is P*N
0060 
0061 % signal length
0062 sig_length = 2^15; % 128*128;
0063 Np = N; % interval of kalman-type signal evolution
0064 
0065 % signal extension for prediction of new coefficients in the streaming window.
0066 eType = 'per';
0067 
0068 % signal dynamics
0069 dType = 'static'; % type of dynamics 'crshift', or 'static'
0070 cshift = -1;
0071 rshift_max = 0.5;
0072 rshift_h = @(z) (rand-0.5)*rshift_max*2;
0073 
0074 % DWT parameters
0075 % type of scaling function
0076 % depth of scaling functions (number of wavelet levels)
0077 % type of extension for DWT (0 - periodic extension, 3 - streaming)
0078 %
0079 %-------------------------------------------------------------------------
0080 % COMPARE recovery using block DWT (sym=0) and overlapping DWT (sym = 3)
0081 %-------------------------------------------------------------------------
0082 wType = 'daub4'; sym = 0;
0083 % wType = 'daub4'; sym = 3;
0084 J = log2(N)-3; 
0085 
0086 % rank-1 update mode
0087 delx_mode = 'mil'; % mil or qr
0088 
0089 % add snapshots of the signal in streaming window and average them before comitting to the output.
0090 avg_output = 0; 
0091  
0092 verbose = 0;
0093 
0094 
0095 %% SM: Sampling modes
0096 % % LM: length of measurement window
0097 % LM = 2*N; % 'universal' sampling scheme (align before the overlapping regions of DWT windows that are measured)
0098 if LM > N
0099     LeftEdge_trunc = 1;
0100 else
0101     LeftEdge_trunc = 0;
0102 end
0103 LeftProj_cancel = 1;
0104 
0105 fprintf('Streaming signal measureed and reconstructed with overlapping windows using DWT ... \n');
0106 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,R) = %d, %d, %d, P = %d, LM = %d, \n wType-%s, J = %d, sym = %d, specified signal-length = %d, \n dType-%s, cshift = %d, rshift_max = %3.4g, Np = %d.', mType, sType, SNR, N, round(N/R), R, P, LM, wType, J, sym, sig_length, dType, cshift, rshift_max, Np);
0107 disp(str0);
0108 
0109 %% DWT setup
0110 % DWT parameters
0111 % Length of each window is L. (extend to adaptive/dyadic windows later?)
0112 % wType = 'daub4'; % type of scaling function
0113 % J = 3; % depth of scaling functions (number of wavelet levels)
0114 % sym = 3; % type of extension for DWT (0 - periodic extension, 3 - streaming)
0115 
0116 in_Psi = []; in_Psi.N = N; in_Psi.J = J; in_Psi.wType = wType; in_Psi.sym = sym;
0117 Psi = create_DWT(in_Psi); % DWT synthesis matrix over a window
0118 L = size(Psi,1);
0119 
0120 %% Signal generation
0121 
0122 % Setup dynamical model
0123 % At every time instance, add to the original/previous signal
0124 % an integer circular shift that is known
0125 % a random drift that is unknown
0126 if strcmpi(dType,'crshift')    
0127     % Generate a random signal
0128     in = []; in.type = sType; in.randgen = 0; in.take_fwt = 0; 
0129     [x_init sig wave_struct] = genSignal(N,in);
0130     
0131     F_h = @(x,cshift,rshift) interp1(1:N,circshift(x,cshift),[1:N]+rshift,'linear','extrap')';
0132     
0133     F0 = zeros(N);
0134     for ii = 1:Np;
0135         F0(:,ii) = F_h(circshift([1; zeros(N-1,1)],ii-1),cshift,0);
0136     end
0137     maxsim = round(sig_length)/N;
0138     sigt = sig; sig = [];
0139     for ii = 1:maxsim;
0140         rshift = rshift_h(1);
0141         sigt = F_h(sigt, cshift, rshift);
0142         sig = [sig; sigt];
0143     end
0144 else
0145     % Generate a random signal
0146     in = []; in.type = sType; in.randgen = 0; in.take_fwt = 0; in.Np = Np;
0147     [x_init sig wave_struct] = genSignal(sig_length,in);
0148 end
0149 sig = [zeros(L-N,1); sig]; 
0150 sig_length = length(sig);
0151  
0152 % view DWT coefficients...
0153 alpha_vec = apply_DWT(sig,N,wType,J,sym);
0154 figure(123); 
0155 subplot(211); imagesc(reshape(alpha_vec,N,length(alpha_vec)/N)); 
0156 axis xy;
0157 subplot(212); plot(alpha_vec); 
0158  
0159 
0160 %% Save results..
0161 streaming_iter = ceil(length(sig)/N);
0162 SIM_stack = cell(streaming_iter,1);
0163 SIM_memory = cell(streaming_iter,1);
0164 
0165 x_vec = zeros(N*streaming_iter,1);
0166 xh_vec = zeros(N*streaming_iter,3);
0167 sig_vec = zeros(length(sig),1);
0168 sigh_vec = zeros(length(sig),3);
0169 
0170 %% Create analysis/synthesis matrix explicitly
0171 in = [];
0172 in.P = P; in.Psi = Psi;
0173 % in.P = P; in.Jp = Jp; in.wType = wType; in.Np = Np; in.sym = sym;
0174 PSI = create_PSI_DWT(in);
0175 
0176 %% Setup sensing matrices
0177 in = []; in.type = mType; in.parts = 4;
0178 genAmat_h = @(M,N) genAmat(M,LM,in);
0179 in.P = P-(LM-N)/N;
0180 in.LM = LM; in.M = M; in.N = N;
0181 PHI = create_PHI(in);
0182 
0183 %% Overall streaming system matrix
0184 T_length = size(PSI,1);
0185 t_ind = 1:T_length;
0186 
0187 sigt = sig(t_ind); % Signal under the LOT window at time t
0188 if sym == 1 || sym == 2
0189     x = pinv(PSI'*PSI)*(PSI'*sigt); % Sparse LOT coefficients
0190 else
0191     x = PSI'*sigt;
0192 end
0193 
0194 PSI_M = PSI(1:end-(L-N),:);
0195 A = PHI*PSI_M;
0196 y = PHI*sigt(1:end-(L-N));
0197 
0198 leny = length(y);
0199 sigma = sqrt(norm(y)^2/10^(SNR/10)/leny);
0200 e = randn(leny,1)*sigma;
0201 y = y+e;
0202 
0203 % parameter selection
0204 % tau = sigma*sqrt(log(N));
0205 tau = max(1e-2*max(abs(A'*y)),sigma*sqrt(log(P*N)));
0206 
0207 maxiter = 2*P*N;
0208 err_fun = @(z) (norm(x-z)/norm(x))^2;
0209 
0210 %% Initialize by solving a rwt L1 problem
0211 in = [];
0212 in.tau = tau; W = tau;
0213 in.W = W;
0214 in.delx_mode = delx_mode;
0215 in.debias = 0;
0216 in.verbose = 0;
0217 in.plots = 0;
0218 in.record = 1;
0219 in.err_fun = err_fun;
0220 tic
0221 for wt_itr = 1:5
0222     
0223     out = l1homotopy(A,y,in);
0224     xh = out.x_out;
0225     iter_bpdn = out.iter;
0226     time_bpdn = toc;
0227     gamma_bpdn = out.gamma;
0228     
0229     % Update weights
0230     xh_old = xh;
0231     
0232     alpha = 1; epsilon = 1;
0233     beta = M*(norm(xh,2)/norm(xh,1))^2;
0234     W = tau/alpha./(beta*abs(xh)+epsilon);
0235     
0236     W_old = W;
0237     yh = A*xh;     
0238     Atr = A'*(A*xh-y);
0239     u =  -W.*sign(xh)-Atr;
0240     pk_old = Atr+u;
0241     
0242     in = out;
0243     in.xh_old = xh;
0244     in.pk_old = pk_old;
0245     in.u = u;
0246     in.W_old = W_old;
0247     in.W = W;
0248 end
0249 W = W_old;
0250 
0251 sim = 1;
0252 x_vec((sim-1)*N+1:sim*N,1) = x(1:N);
0253 xh_vec((sim-1)*N+1:sim*N,1:3) = [xh(1:N) xh(1:N) xh(1:N)]; 
0254 
0255 s_ind = 1:L;
0256 sig_temp = Psi*xh(1:N);
0257 s_offset = 0;
0258 
0259 sig_vec(s_ind) = sigt(s_ind);
0260 sigh_vec(s_ind,1:3) = sigh_vec(s_ind,1:3)+[sig_temp sig_temp sig_temp];
0261 
0262 % average instantaneous estimates before committing to output...
0263 estimate_buffer = repmat(xh(1:(P-1)*N,1),1,P-1)/(P-1);
0264 
0265 xh_streamingRWT = xh;
0266 x_sparsa = xh;
0267 x_yall1 = xh;
0268 
0269 done = 0;
0270 while ~done
0271      
0272     %% Update the solution after updating the measurement matrix and/or the
0273     % sparse signal
0274     y_old = y; x_old = x;
0275     
0276     sigt_old = sigt; t_ind_old = t_ind;
0277     PHI_old = PHI; 
0278     
0279     % Shift the sampling window
0280     t_ind = t_ind+N;
0281     if t_ind(end) > length(sig)
0282         break;
0283     end
0284     sigt = sig(t_ind);
0285     
0286     % System matrix setup...
0287     % Shift up and left
0288     PHI(1:end-M,1:end-N) = PHI(M+1:end,N+1:end);
0289     % new measurement matrix
0290     Phi = genAmat_h(M,LM);
0291     PHI(end-M+1:end,end-LM+1:end) = Phi;
0292     
0293     
0294     A = PHI*PSI_M;
0295     y = PHI*sigt(1:end-(L-N));
0296     
0297     % Sparse LOT coeffs. irrespective of SM
0298     if sym == 1 || sym == 2
0299         x = pinv(PSI'*PSI)*(PSI'*sigt);
0300     else
0301         x = PSI'*sigt;
0302     end
0303     
0304     % shift old measurements and add one new set of measurementts
0305     e(1:end-M) = e(M+1:end);
0306     e(end-M+1:end) = randn(M,1)*sigma;
0307     y= y+e;
0308     
0309     A0 = A; x0 = x; y0 = y;
0310     for solver = {'l1homotopy','sparsa','yall1'}
0311         solver = char(solver);
0312         switch solver
0313             case 'l1homotopy'
0314                 xh = xh_streamingRWT;
0315             case 'sparsa'
0316                 xh = x_sparsa;
0317             case 'yall1'
0318                 xh = x_yall1;
0319         end
0320         y = y0; A = A0; x = x0;
0321         xh_old = xh;
0322         
0323           
0324         % REMOVE the part of outgoing DWT projection in the overlapping region
0325         % on left side of streaming window...
0326         if LeftProj_cancel
0327             y(1:M) = y(1:M)-PHI(1:M,1:(L-N))*(Psi(end-(L-N)+1:end,:)*xh_old(1:N));
0328         end
0329         
0330         %% Update the signal estimate (for warm start)
0331         xh(1:end-N) = xh(N+1:end);
0332         % truncate small values in xh
0333         % xh(abs(xh)<tau/sqrt(log(P*N))) = 0;
0334         
0335         % --------------------------------------------%
0336         %  Linear prediction for the incoming samples %
0337         % --------------------------------------------%
0338         %
0339         %         xh1 = xh(end-N+1:end);
0340         %
0341         %         if strcmpi(eType,'per')
0342         %             % periodic extension
0343         %             xh1 = xh(end-N+1:end);
0344         %         else
0345         %             sig_temp = PSI_M(1:end-N,1:(P-1)*N)*xh(1:end-N);
0346         %             sig_temp(1:L-N) = sig_temp(1:L-N)+Psi(end-(L-N)+1:end,:)*xh_old(1:N);
0347         %             % if sim == 1
0348         %             %     fprintf('Extension type: %s \n',eType);
0349         %             % end
0350         %             switch eType
0351         %                 case 'asym'
0352         %                     % anti-symmetric extension
0353         %                     sig_temp2 = [sig_temp; -sig_temp(end:-1:1); sig_temp];
0354         %                 case 'sym'
0355         %                     % symmetrci extension
0356         %                     sig_temp2 = [sig_temp; sig_temp(end:-1:1); sig_temp];
0357         %                 otherwise
0358         %                     disp('cant do it sire');
0359         %             end
0360         %             sig_temp2 = sig_temp2(1:size(PSI,1));
0361         %             if sym == 1 || sym == 2
0362         %                 xh_temp = pinv(PSI'*PSI)*(PSI'*sig_temp2);
0363         %             else
0364         %                 xh_temp = PSI'*sig_temp2;
0365         %             end
0366         %             xh1 = xh_temp(end-N+1:end);
0367         %         end
0368         %
0369         %         gamma1 = find(abs(xh1)>tau);
0370         %         [val ind] = sort(abs(xh1),'descend');
0371         %         gamma1 = ind(1:min(length(gamma1),ceil(M/log2(N))));
0372         %
0373         %         yt = y(end-M+1:end)-A(end-M+1:end,1:end-N)*xh(1:end-N);
0374         %         At = A(end-M+1:end,end-N+1:end);
0375         %         gamma3 = gamma1;
0376         %         %
0377         %         % NOTE TODO: Need to replace prediction with rank-update...
0378         %         %
0379         %         % gamma_all = [find(xh(1:end-N)); gamma3+(P-1)*N];
0380         %         % xh_all = 0*xh;
0381         %         % xh_all(gamma_all) = A(:,gamma_all)\y; % pinv(A(:,gamma_all)'*A(:,gamma_all))*(A(:,gamma_all)'*y-W(gamma_all).*sign(xh(gamma_all)));
0382         %
0383         %         % if SM == 3
0384         %         %     Aty = At'*yt;
0385         %         %     [val ind] = sort(abs(Aty),'descend');
0386         %         %     gamma2 = ind(1:min(length(gamma1),ceil(M/4)));
0387         %         %     gamma3 = union(gamma1,gamma2);
0388         %         % end
0389         %         %
0390         %         xh_t1 = 0*xh(end-N+1:end);
0391         %         xh_t1(gamma3) = pinv(At(:,gamma3)'*At(:,gamma3)+2*tau*eye(length(gamma3)))*(At(:,gamma3)'*yt);
0392         %         xh_t1(abs(xh_t1)<tau/sqrt(log(P*N))) = 0;
0393         %         xh(end-N+1:end) = xh_t1;
0394         
0395         
0396         % Truncate small coefficients... ?
0397         xh(abs(xh)<tau/sqrt(log(P*N))) = 0;
0398         fig(111);
0399         plot([x xh]);
0400         
0401         % Remove the top-left edge of the system matrix
0402         if LeftEdge_trunc
0403             % fprintf('Consider oldest set of LOT coefficients to be fully known, and remove their contribution from the measurements... \n');
0404             alpha0h = xh(1:N);
0405             xh = xh(N+1:end);
0406             y = y-A(:,1:N)*alpha0h;
0407             
0408             A = A(:,N+1:end);
0409             alpha0 = x(1:N);
0410             x = x(N+1:end);
0411         end
0412         
0413         % Update weights
0414         epsilon = 1;
0415         beta = max(epsilon,M*(norm(xh,2)/norm(xh,1))^1);
0416         W = tau/alpha./(beta*abs(xh)+epsilon);
0417         
0418         W_old = W;
0419         
0420         
0421         if strcmpi(solver,'l1homotopy')
0422             
0423             homotopy_mode = 'dummy';
0424             switch homotopy_mode
0425                 case 'dummy'
0426                     % create a dummy variable...
0427                     % use homotopy on the measurements...
0428                     % in principle, we can start with any xh with this formulation
0429                     % and any starting value of tau or W...
0430                     gamma = find(xh);
0431                     M_trunc = size(A,1); % P*(M-1);
0432                     % M_trunc = round(P*(N-1)/5);
0433                     if length(gamma) >= M_trunc
0434                         disp('length of gamma exceeded number of rows');
0435                         [xh_sort ind_sort] = sort(abs(xh),'descend');
0436                         xh(ind_sort(M_trunc+1:end)) = 0;
0437                         gamma = ind_sort(1:M_trunc);
0438                     end
0439                     Atr = A'*(A*xh-y);
0440                     u =  -W.*sign(xh)-Atr;
0441                     pk_old = Atr+u;
0442                 otherwise
0443                     didp('Go back ... no escape');
0444             end
0445             
0446             
0447             in = out;
0448             gamma_old = gamma;
0449             in.gamma = gamma_old;
0450             switch delx_mode
0451                 case 'mil';
0452                     % in.delx_mode = 'mil';
0453                     % The following gram matrix and its inverse can be used from the
0454                     % previous homotopy. Too lazy to include that right now...
0455                     % wt BPDN homotopy update
0456                     AtAgx = A(:,gamma_old)'*A(:,gamma_old);
0457                     iAtAgx = pinv(AtAgx);
0458                     in.iAtA = iAtAgx;
0459                 case {'qr','chol'};
0460                     % in.delx_mode = 'qr';
0461                     [Q R] = qr(A(:,gamma_old),0);
0462                     in.Q = Q; in.R = R;
0463                 case 'qrM'
0464                     % in.delx_mode = 'qrM';
0465                     [Q0 R0] = qr(A(:,gamma_old));
0466                     in.Q0 = Q0; in.R0 = R0;
0467             end
0468             
0469             in.xh_old = xh;
0470             in.pk_old = pk_old;
0471             in.u = u;
0472             in.W = W;
0473             in.delx_mode = delx_mode;
0474             in.debias = 0;
0475             in.verbose = 0;
0476             in.plots = 0;
0477             in.record = 1;
0478             in.err_fun = @(z) (norm(x-z)/norm(x))^2;
0479             tic
0480             out = l1homotopy(A,y,in);
0481             time_streamingRWT = toc;
0482             xh_streamingRWT = out.x_out;
0483             gamma_streamingRWT = out.gamma;
0484             iter_streamingRWT = out.iter;
0485             % Reconstructed signal
0486             if LeftEdge_trunc
0487                 x = [alpha0; x];
0488                 xh_streamingRWT = [alpha0h; xh_streamingRWT];
0489             end         
0490         elseif  strcmpi(solver,'sparsa')
0491             %% SpaRSA
0492             x_sparsa = xh; W_sparsa = W/tau; iter_sparsa = 0; time_sparsa = 0;
0493             if norm(y) > 1e-3
0494                 psi_function = @(x,tau) soft(x,tau*W_sparsa);
0495                 phi_function = @(x) sum(abs(W_sparsa.*x));
0496                 tic;
0497                 [x_sparsa,x_debias_SpaRSA_m,obj_SpaRSA_m_cont,...
0498                     times_SpaRSA_m_cont,debias_start_SpaRSA_m,mse_SpaRSA_m,taus_SpaRSA_m, numA, numAt]= ...
0499                     SpaRSA_adpW(y,A,tau,...
0500                     'Monotone',0,...
0501                     'adp_wt',0,...
0502                     'W_new',W_sparsa,...
0503                     'Debias',0,...
0504                     'Initialization',x_sparsa,...
0505                     'StopCriterion',2,...
0506                     'ToleranceA',1e-4,...
0507                     'psi',psi_function,...
0508                     'phi',phi_function,...
0509                     'Safeguard',1,...
0510                     'MaxiterA',maxiter,...
0511                     'Verbose',0,...
0512                     'True_x',x,...
0513                     'Continuation',1,...
0514                     'Continuationsteps',-1);
0515                 
0516                 time_sparsa = toc;
0517                 iter_sparsa = (numA+numAt)/2;
0518                 error_sparsa = norm(x-x_sparsa)/norm(x);            
0519             end
0520             % Reconstructed signal
0521             if LeftEdge_trunc
0522                 x = [alpha0; x];
0523                 x_sparsa = [alpha0h; x_sparsa];
0524             end
0525         elseif strcmpi(solver,'yall1')            
0526             
0527             %% YALL1
0528             % set options
0529             digit = 4; if sigma > 0; digit = 4; end
0530             opts = [];
0531             opts.tol = 10^(-digit);
0532             opts.weights = W/tau;
0533             opts.print = 0;
0534             opts.maxit = maxiter;
0535             opts.nonorth = 1;
0536             % opts.x0 = xh;
0537             opts.nu = 0; opts.rho = tau;
0538             tic;
0539             [x_yall1,Out_yall1] = yall1(A,y,opts);
0540             % time_yall1 = [time_yall1 Out.cputime];
0541             time_yall1 = toc;
0542             iter_yall1 = (Out_yall1.cntA+Out_yall1.cntAt)/2;
0543             err_yall1 = norm(x-x_yall1)/norm(x);
0544             % Reconstructed signal
0545             if LeftEdge_trunc
0546                 x = [alpha0; x];
0547                 x_yall1 = [alpha0h; x_yall1];
0548             end         
0549             
0550             if max(abs(x_yall1)) > 500
0551                 stp = 1;
0552             end
0553         end 
0554     end
0555     %% Record results
0556     sim = sim+1;
0557     SIM_stack{sim} = [sim, tau, ...
0558         norm(x-xh_streamingRWT)^2/norm(x)^2, sum(iter_streamingRWT,2), sum(time_streamingRWT,2), ...
0559         norm(x-x_sparsa)^2/norm(x)^2, sum(iter_sparsa,2), sum(time_sparsa,2), ...
0560         norm(x-x_yall1)^2/norm(x)^2, sum(iter_yall1,2), sum(time_yall1,2)];
0561     
0562     % print and plot
0563     if mod(sim-1,verbose) == 0 && verbose
0564         fprintf('streaming iter. %d. tau = %3.4g, (err,iter,time): streamingRWT homotopy-%3.4g,%3.4g,%3.4g, SpaRSA-%3.4g,%3.4g,%3.4g, YALL1-%3.4g,%3.4g,%3.4g. \n', ...
0565             SIM_stack{sim});
0566     end
0567     
0568     %% Plot DWT coeffs. on the window
0569     fig(1); subplot(211);
0570     plot([x xh_streamingRWT x_sparsa x_yall1]);
0571     title('Comparison betweeen the original and reconstructed signal')
0572     
0573     %% Reconstructed signal
0574     x_vec((sim-1)*N+1:sim*N,1) = x(1:N);
0575     xh = xh_streamingRWT;
0576     
0577     % remove the oldest estimate, shift the remaining up and left, and add the new estimate
0578     estimate_buffer = [[estimate_buffer(N+1:end,2:end); zeros(N,P-2)] xh(1:end-N)/(P-1)];
0579     if avg_output 
0580         xh_est = xh(1:N);
0581         xh(1:N) = sum(estimate_buffer(1:N,:),2);
0582         % fig(123); plot([xh_est xh(1:N) x(1:N)])
0583         if sim == 2
0584             disp('output is averaged');
0585         end
0586     end
0587     xh_vec((sim-1)*N+1:sim*N,1) = xh(1:N);
0588     xh_vec((sim-1)*N+1:sim*N,2) = x_sparsa(1:N);
0589     xh_vec((sim-1)*N+1:sim*N,3) = x_yall1(1:N);
0590     
0591     
0592     s_ind = t_ind(1:L);
0593     sig_vec(s_ind) = sigt(1:L);
0594     sigh_vec(s_ind,1) = sigh_vec(s_ind,1)+Psi*xh(1:N);
0595     sigh_vec(s_ind,2) = sigh_vec(s_ind,2)+Psi*x_sparsa(1:N);
0596     sigh_vec(s_ind,3) = sigh_vec(s_ind,3)+Psi*x_yall1(1:N);
0597 
0598     
0599     
0600     %% plot recovered signals
0601     fig(1); subplot(212)
0602     plot([sig_vec(1:s_ind(end)) sigh_vec(1:s_ind(end),1)]);
0603  
0604     drawnow;
0605 end
0606 
0607 mS =  sum(cell2mat(SIM_stack),1);
0608 fprintf('Summed results: streaming_iter %d. tau = %3.4g, \n solver-(err,iter,time): \n streamingRWT homotopy-%3.4g,%3.4g,%3.4g; \n SpaRSA-%3.4g,%3.4g,%3.4g; \n YALL1-%3.4g,%3.4g,%3.4g. \n', streaming_iter, mS(2:end));
0609 mS =  mean(cell2mat(SIM_stack),1);
0610 fprintf('Average results: streaming_iter %d. tau = %3.4g, \n solver-(err,iter,time): \n streamingRWT homotopy-%3.4g,%3.4g,%3.4g; \n SpaRSA-%3.4g,%3.4g,%3.4g; \n YALL1-%3.4g,%3.4g,%3.4g. \n', streaming_iter, mS(2:end));
0611 
0612 % l1homotopy-%3.4g,%3.4g,%3.4g;
0613 st_ind = 2*N;
0614 err_l1homotopy = norm(sig_vec(st_ind+1:N*sim)-sigh_vec(st_ind+1:N*sim,1))^2/norm(sig_vec(st_ind+1:N*sim))^2;
0615 err_sparsa = norm(sig_vec(st_ind+1:N*sim)-sigh_vec(st_ind+1:N*sim,2))^2/norm(sig_vec(st_ind+1:N*sim))^2;
0616 err_yall1 = norm(sig_vec(st_ind+1:N*sim)-sigh_vec(st_ind+1:N*sim,3))^2/norm(sig_vec(st_ind+1:N*sim))^2;
0617 % err_l1homotopy = norm(sig_vec(L-N+1:N*sim)-sigh_vec(L-N+1:N*sim,1))^2/norm(sig_vec(L-N+1:N*sim))^2;
0618 % err_sparsa = norm(sig_vec(L-N+1:N*sim)-sigh_vec(L-N+1:N*sim,2))^2/norm(sig_vec(L-N+1:N*sim))^2;
0619 % err_yall1 = norm(sig_vec(L-N+1:N*sim)-sigh_vec(L-N+1:N*sim,3))^2/norm(sig_vec(L-N+1:N*sim))^2;
0620 % fprintf('Signal error: l1homotopy-%3.4g, sparsa-%3.4g, yall1-%3.4g.\n',err_l1homotopy,err_sparsa,err_yall1);
0621 fprintf('Signal MSE: l1homotopy-%3.4g, sparsa-%3.4g, yall1-%3.4g.\n',([err_l1homotopy,err_sparsa,err_yall1]));
0622 fprintf('Signal SER (in dB): l1homotopy-%3.4g, sparsa-%3.4g, yall1-%3.4g.\n',-10*log10([err_l1homotopy,err_sparsa,err_yall1]));
0623 
0624 
0625 
0626 %% plot signal and reconstruction error
0627 x_len = min(length(x_vec),length(xh_vec))-(P-1)*N;
0628 sig_len = min(length(sig_vec),length(sigh_vec))-(P-1)*N-L-N;
0629 x_vec1 = x_vec(1:x_len);
0630 xh_vec1 = xh_vec(1:x_len,1);
0631 sig_vec1 = sig_vec(1:sig_len);
0632 sigh_vec1 = sigh_vec(1:sig_len,1);
0633 
0634 fig(123); 
0635 subplot(221); 
0636 plot((1:sig_len)/N,sig_vec1, 'LineWidth',1);
0637 axis tight;
0638 title('original signal')
0639 subplot(2,2,2)
0640 plot((1:sig_len)/N,sigh_vec1-sig_vec1, 'LineWidth',1);
0641 axis tight
0642 title('reconstruction error')
0643 subplot(2,2,3); 
0644 imagesc(reshape(x_vec1,N,x_len/N)); axis xy;
0645 axis tight;
0646 title('DWT coefficients');
0647 colorbar 
0648 subplot(2,2,4); 
0649 imagesc(reshape(x_vec1-xh_vec1,N,x_len/N),[0 max(abs(x_vec1))/20]); axis xy
0650 axis tight
0651 title('reconstruction error (DWT)');
0652 colorbar
0653 
0654         
0655 fig(3);
0656 % view DWT coefficients...
0657 alpha_vec1 = apply_DWT(sig_vec1,N,wType,J,sym);
0658 alphah_vec1 = apply_DWT(sigh_vec1,N,wType,J,sym);
0659 subplot(221); imagesc(reshape(alpha_vec1,N,length(alpha_vec1)/N));
0660 axis xy;
0661 title('original')
0662 subplot(222); imagesc(reshape(alphah_vec1,N,length(alpha_vec1)/N));
0663 axis xy;
0664 title('reconstructed');
0665 subplot(212); plot([alpha_vec1 alphah_vec1]); 
0666 title('comparison');