Compressed sensing STA low memory version

cssta_lowmem.m

%Download and uncompress:
%http://users.ece.gatech.edu/~justin/l1magic/downloads/l1magic-1.11.zip

addpath l1magic/

%W with smooth connected points
L = 32;
z = linspace(-.25,1.25,L);
[Z1,Z2] = meshgrid(z,z);
rho = .015;
R = sqrt(Z1.*Z1+Z2.*Z2);
W = (1-rho)<R & R <(1+rho);
W = 10*W.*(1+Z1+Z2);

%set up params
T = 2^16;
thresh = 10;
q = .01;%percent flips

%set up random projection
N = length(W(:));
K = 256;
A = rand(K,N);
A = orth(A')';

%set up vectors
spikes = zeros(1,T);
W_sta = zeros(size(W));
W_comp_sta = zeros(K,1);

%run all
for t = 1:T
    stim = randn([L,L]);
    r = sum(sum(stim.*W));
    spikes(t) = 2*(thresh<r)-1;
    spikes(t) = .5*(spikes(t)*(2*(rand()>q)-1)+1);
    
    %accum STA estimate
    W_sta= W_sta + spikes(t)*stim;
    
    %project stimulus and accum compressed STA
    unrolled_stim_proj = A*stim(:);
    W_comp_sta =W_comp_sta+spikes(t)*unrolled_stim_proj;
end

%normalize sta estimates
W_sta = W_sta/sum(spikes);
W_comp_sta = W_comp_sta/sum(spikes);



%Uncompress STA estimate
x0 = A'*W_comp_sta ;
xp = l1eq_pd(x0, A, [], W_comp_sta, 2e-3);
W_cssta = reshape(xp,size(W));

%results
figure(1)
subplot(131)
imagesc(W)
title('Exact filter')
axis image
subplot(132)
imagesc(W_sta)
title('STA estimated filter')
axis image
subplot(133)
imagesc(W_cssta)
title('Uncompressed CCSTA estimated filter')
axis image