demo_streamingLOT

 demo_streamingLOT

 Solves the following streaming problem
 min_x  \|x\|_1 + 1/2*||A  x - y||_2^2

 by adding and removing one set of measurements at every time instant

 Applications:
           streaming signal recovery using lapped orthogonal transform (LOT)

 We can also add any other regularization operator in the reconstruction

 Written by: Salman Asif, Georgia Tech
 Email: sasif@gatech.edu
 Created: August 2012

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