Quick demo of compressed sensing Spike triggered average analysis

csstademo.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^19;
stim = randn([L,L,T]);
thresh = 10;
q = .01;%percent flips

%run stimulus response experiment 
spikes = zeros(1,T);
for t = 1:T
    r = sum(sum(stim(:,:,t).*W));
    spikes(t) = 2*(thresh<r)-1;
end
spikes = spikes.*(2*(rand(size(spikes))>q)-1);
spikes = .5*(spikes+1);

%estimate STA using full stimulus
W_sta = zeros(size(W));
for t=1:T
    W_sta= W_sta + spikes(t)*stim(:,:,t);
end
W_sta = W_sta/sum(spikes);

%Now repeat STA in projected stimulus space

%Construct projections
N = length(W(:));
K = 256;
A = rand(K,N);
A = orth(A')';

%Projected stimulus observations
unrolled_stim_proj = zeros(K,T);
for t=1:T
    stim_t = stim(:,:,t);
    unrolled_stim_proj(:,t) = A*stim_t(:);
end

%estimate STA using compressed stimulus
W_comp_sta = zeros(K,1);
for t=1:T
    W_comp_sta =W_comp_sta+spikes(t)*unrolled_stim_proj(:,t);
end
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 CSTA estimated filter')
axis image