03-27-13

block_levinson.m

@@ -0,0 +1,41 @@
+function x = block_levinson(y, L)
+%BLOCK_LEVINSON Block Levinson recursion for efficiently solving 
+%symmetric block Toeplitz matrix equations.
+%   BLOCK_LEVINSON(Y, L) solves the matrix equation T * x = y, where T 
+%   is a symmetric matrix with block Toeplitz structure, and returns the 
+%   solution vector x. The matrix T is never stored in full (because it 
+%   is large and mostly redundant), so the input parameter L is actually 
+%   the leftmost "block column" of T (the leftmost d columns where d is 
+%   the block dimension).
+
+%   Author: Keenan Pepper
+%   Last modified: 2007-12-23
+
+%   References:
+%     [1] Akaike, Hirotugu (1973). "Block Toeplitz Matrix Inversion".
+%     SIAM J. Appl. Math. 24 (2): 234-241
+
+s = size(L);
+d = s(2);                 % Block dimension
+N = s(1) / d;             % Number of blocks
+
+B = reshape(L, [d,N,d]);  % This is just to get the bottom block row B
+B = permute(B, [1,3,2]);  % from the left block column L
+B = flipdim(B, 3);
+B = reshape(B, [d,N*d]);
+
+f = L(1:d,:)^-1;          % "Forward" block vector
+b = f;                    % "Backward" block vector
+x = f * y(1:d);           % Solution vector
+
+for n = 2:N
+    ef = B(:,(N-n)*d+1:N*d) * [f;zeros(d)];
+    eb = L(1:n*d,:)' * [zeros(d);b];
+    ex = B(:,(N-n)*d+1:N*d) * [x;zeros(d,1)];
+    A = [eye(d),eb;ef,eye(d)]^-1;
+    fn = [[f;zeros(d)],[zeros(d);b]] * A(:,1:d);
+    bn = [[f;zeros(d)],[zeros(d);b]] * A(:,d+1:end);
+    f = fn;
+    b = bn;
+    x = [x;zeros(d,1)] + b * (y((n-1)*d+1:n*d) - ex);
+end

hippo/cfastica_defl.m

@@ -0,0 +1,28 @@
+function W = cfastica_defl(x)
+eps = 0.1; % epsilon in G
+[x,whiteningMatrix,dewhiteningMatrix,zerophaseMatrix] = whiten(x,.1);
+n = size(x,1);
+
+W = zeros(n);
+maxcounter = 40;
+for k = 1:n
+    w = rand(n,1) + 1i*rand(n,1);
+    counter = 0;
+    wold = zeros(n,1);
+    while min(sum(abs(abs(wold) - abs(w))), maxcounter - counter) > 0.001;
+        wold = w;
+        g = 1./(eps + abs(w'*x).^2);
+        dg = -1./(eps + abs(w'*x).^2).^2;
+        w = mean(x .* (ones(n,1)*conj(w'*x)) .* (ones(n,1)*g), 2) - ...
+            mean(g + abs(w'*x).^2 .* dg) * w;
+        w = w / norm(w);
+        % Decorrelation:
+        w = w - W*W'*w;
+        w = w / norm(w);
+        counter = counter + 1;
+    end
+    W(:,k) = w;
+end
+W = W';
+x = W * x;
+figure;imagesc(abs(x));

hippo/cfastica_defl.m~

@@ -0,0 +1,30 @@
+function W = cfastica_defl(x)
+eps = 0.1; % epsilon in G
+[x,whiteningMatrix,dewhiteningMatrix,zerophaseMatrix] = whiten(mixedsig);
+n = size(x,1);
+
+  W = zeros(n,n);
+  maxcounter = 40;
+  for k = 1:n
+    w = rand(n,1) + i*rand(n,1);
+    counter = 0;
+    wold = zeros(n,1);
+    while min(sum(abs(abs(wold) - abs(w))), maxcounter - counter) > 0.001;
+      wold = w;
+      g = 1./(eps + abs(w'*x).^2);
+      dg = -1./(eps + abs(w'*x).^2).^2;
+      w = mean(x .* (ones(n,1)*conj(w'*x)) .* (ones(n,1)*g), 2) - ...
+          mean(g + abs(w'*x).^2 .* dg) * w;
+      w = w / norm(w);
+      % Decorrelation:
+      w = w - W*W'*w;
+      w = w / norm(w);
+      counter = counter + 1;
+      G = log(eps + abs(w'*x).^2);
+    end
+    W(:,k) = w;
+    counters(k) = counter;
+  end; 
+
+W_hat = W';
+e_hat = W_hat * x;

hippo/findPlaces.m

@@ -0,0 +1,71 @@
+function [fits cis acts] = findPlaces(fields,thresh,sign)
+
+buffer = floor(size(fields,3)/2/10);
+if ~exist('sign','var') || sign > 0
+    sign = 1;
+end
+%[xs ys] = meshgrid(1:size(fields,2),1:size(fields,3));
+for i = 1:size(fields,1)
+    fits{i} = [];cis{i} = [];
+    for j = 1:2
+        tempFull = sign*squeeze(fields(i,:,(j-1)*size(fields,3)/2+(1:size(fields,3)/2)));
+        %tempFull = filtfilt(gausswin(buffer/2),sum(gausswin(buffer/2)),tempFull')';
+        temp = mean(tempFull);
+        tempF = mean(abs(tempFull));
+        [~,locs] = findpeaks(mean(tempF,1),'minpeakheight',thresh,'minpeakdistance',buffer);
+        x0 = [];lb = [];ub = [];
+        for k = 1:numel(locs)
+            x0 = [x0; 1/20 -.1 locs(k) real(temp(locs(k)))' imag(temp(locs(k)))'];
+            lb = [lb; 0 -inf locs(k)-buffer/2 -inf -inf];
+            ub = [ub; inf 0 locs(k)+buffer/2 inf inf];
+            %allMags(k,:) = tempF(:,locs(k))'/norm;
+        end
+        subplot(311);
+        imagesc(complexIm(tempFull,[],[],[],[],thresh*2));title((2*(i-1)+j));
+        hold all;scatter(locs,ones(numel(locs),1)*size(tempFull,2)/2,'w','filled');hold off;
+        if numel(locs)
+            %f = @(x,xdata) myfun(x,xdata,allMags);
+            [beta,~,resid,~,~,~,J] = lsqcurvefit(@myfun,x0,size(temp,2),[real(temp);imag(temp)],lb,ub);
+            ci = nlparci(beta,resid,'jacobian',J);
+            ci = reshape(ci(:,2)-ci(:,1),size(beta));
+            for k = 1:numel(locs)
+                mor = makeCmor(size(temp,2),beta(k,1),beta(k,2),beta(k,3));
+                acts{i}(k,:) = (mor*tempFull')/(mor*mor');
+            end
+            c = myfun(beta,size(temp,2),acts{i});%,allMags
+            c = complex(c(1:size(c,1)/2,:),c(size(c,1)/2+1:end,:));
+            subplot(312);imagesc(complexIm(c,[],[],[],[],thresh*2));hold all;
+            scatter(beta(:,3),ones(size(beta,1),1)*size(tempFull,2)/2,'w','filled');hold off;title(numel(locs));drawnow;
+            beta(:,3) = beta(:,3) + (j-1)*size(fields,3)/2;
+            fits{i} = [fits{i} beta'];
+            cis{i} = [cis{i} ci'];
+            %fits(2*(i-1)+j).x
+            tc = tempFull-c;
+            tt = [mean(tempFull);mean(c)];
+            subplot(313);sPlot(complex(abs(tt)/max(abs(tt(:)))*6,angle(tt)),[],0);%imagesc(complexIm(tc,[],[],[],[],thresh*2));
+            drawnow;input('');%pause(1);
+            %[sum(mean(tempFull.*conj(tempFull))) sum(mean(c.*conj(c))) sum(mean(tc.*conj(tc)))]
+        end
+        %ac = 0;
+        %for k = 1:size(temp,1)
+        %    ac = ac + xcov(temp(k,:),buffer);
+        %end
+        %
+    end
+end
+
+function F = myfun(x,len,mags)%
+F = 0;
+for i = 1:size(x,1);
+    thisMor = makeCmor(len,x(i,1),x(i,2),x(i,3));
+    if exist('mags','var')
+        F = F + mags(i,:)'*thisMor;%
+    else
+        F = F + complex(x(i,4),x(i,5))*thisMor;%
+    end
+end
+F = [real(F);imag(F)];
+
+function [psi,X] = makeCmor(len,Fb,Fc,off)
+X = (1:len)-off;  % wavelet support.
+psi = exp(2*1i*pi*Fc*X).*exp(-(X.*X)*Fb);%((pi*Fb)^(-0.5))*

hippo/findPlaces.m~

@@ -0,0 +1,65 @@
+function [fits cis acts] = findPlaces(fields,thresh,sign)
+
+buffer = floor(size(fields,3)/2/10);
+if ~exist('sign','var') || sign > 0
+    sign = 1;
+end
+%[xs ys] = meshgrid(1:size(fields,2),1:size(fields,3));
+for i = 1:size(fields,1)
+    fits{i} = [];cis{i} = [];
+    for j = 1:2
+        tempFull = sign*squeeze(fields(i,:,(j-1)*size(fields,3)/2+(1:size(fields,3)/2)));
+        tempFull = filtfilt(gausswin(buffer),sum(gausswin(buffer)),tempFull')';
+        temp = mean(tempFull);
+        tempF = abs(temp);
+        [~,locs] = findpeaks(mean(tempF,1),'minpeakheight',thresh,'minpeakdistance',buffer);
+        x0 = [];
+        for k = 1:numel(locs)
+            x0 = [x0; 1/20 -.1 locs(k) tempF(locs(k))'];
+            %allMags(k,:) = tempF(:,locs(k))'/norm;
+        end
+        subplot(311);
+        imagesc(complexIm(tempFull,[],[],[],[],thresh*2));title((2*(i-1)+j));
+        if numel(locs)
+            %f = @(x,xdata) myfun(x,xdata,allMags);
+            [beta,~,resid,~,~,~,J] = lsqcurvefit(@myfun,x0,size(temp,2),[real(temp);imag(temp)]);
+            ci = nlparci(beta,resid,'jacobian',J);
+            ci = reshape(ci(:,2)-ci(:,1),size(beta));
+            for k = 1:numel(locs)
+                mor = makeCmor(size(temp,2),beta(k,1),beta(k,2),beta(k,3));
+                acts{i}(k,:) = (mor*tempFull')/(mor*mor');
+            end
+            c = myfun(beta,size(temp,2),acts{i});%,allMags
+            c = complex(c(1:size(c,1)/2,:),c(size(c,1)/2+1:end,:));
+            beta(:,3) = beta(:,3) + (j-1)*size(fields,3)/2;
+            fits{i} = [fits{i} beta'];
+            cis{i} = [cis{i} ci'];
+            %fits(2*(i-1)+j).x
+            subplot(312);
+            imagesc(complexIm(c,[],[],[],[],thresh*2));title(numel(locs));
+            subplot(313);imagesc(complexIm(tempFull-c,[],[],[],[],thresh*2));
+            title(drawnow;
+        end
+        %ac = 0;
+        %for k = 1:size(temp,1)
+        %    ac = ac + xcov(temp(k,:),buffer);
+        %end
+        %
+    end
+end
+
+function F = myfun(x,len,mags)%
+F = 0;
+for i = 1:size(x,1);
+    thisMor = makeCmor(len,x(i,1),x(i,2),x(i,3));
+    if exist('mags','var')
+        F = F + mags(i,:)'*thisMor;%
+    else
+        F = F + x(i,4)*thisMor;
+    end
+end
+F = [real(F);imag(F)];
+
+function [psi,X] = makeCmor(len,Fb,Fc,off)
+X = (1:len)-off;  % wavelet support.
+psi = exp(2*1i*pi*Fc*X).*exp(-(X.*X)*Fb);%((pi*Fb)^(-0.5))*

hippo/findPlacesNon.m

@@ -0,0 +1,68 @@
+function [filts acts] = findPlacesNon(fields,thresh)
+warning off;
+fields(:) = zscore(fields(:));
+buffer = 24;%floor(size(fields,3)/2/10);
+for i = 1:size(fields,1)
+        filts{i} = [];acts{i} = [];
+        tempFull = squeeze(fields(i,:,:));
+        tempFull = [zeros(size(tempFull,1),buffer/2) tempFull zeros(size(tempFull,1),buffer/2)];
+        tempFullF = filtfilt(gausswin(buffer/2),sum(gausswin(buffer/2)),tempFull')';
+        tempF = mean(abs(tempFullF));
+        [~,locs] = findpeaks(tempF,'minpeakheight',thresh,'minpeakdistance',buffer/2);
+        subplot(221);hold off;imagesc(complexIm(tempFull));hold all;scatter(locs,ones(numel(locs),1)*size(tempFull,2)/2,'w','filled');drawnow;
+        if numel(locs)
+            wein = zeros(numel(locs),buffer+1);
+            tempM = mean(tempFull);
+            for k = 1:numel(locs)
+                wein(k,:) = tempM(locs(k)+(-buffer/2:buffer/2));
+            end
+            for m = 1:3
+                if min(size(wein)) == 1
+                    wein = reshape(wein,[ buffer+1 numel(locs)]).';
+                    for k = 1:size(wein,1)
+                        mx = round(sum((1:size(wein,2)).*abs(wein(k,:)))/sum(abs(wein(k,:))));%max(abs(wein(k,:)));
+                        wein(k,:) = circshift(wein(k,:),[0 buffer/2+1-mx]);
+                    end
+                end
+                wein = bsxfun(@rdivide,wein,sqrt(sum(abs(wein).^2,2)));
+                weinOld = wein;
+                toepReg = zeros(numel(locs)*(buffer+1),numel(tempFull));
+                for k = 1:numel(locs)
+                    multReg = zeros(size(tempFull,1),size(tempFull,2));
+                    inds = locs(k)+(-buffer/4:buffer/4);
+                    mxInds = zeros(size(tempFull,1),1);absInds = mxInds;compInds = mxInds;
+                    for l = 1:size(tempFull,1)
+                        tmp = circshift(conv(conj(fliplr(wein(k,:))),tempFull(l,:),'full'),[0 -buffer/2]);
+                        tmp(~ismember(1:numel(tmp),inds)) = 0;
+                        [absInds(l),mxInds(l)] = max(abs(tmp));
+                        if absInds(l) > thresh && mxInds(l) ~= min(inds) && mxInds(l) ~= max(inds)
+                            multReg(l,max(1,mxInds(l))) = tmp(mxInds(l));
+                            compInds(l) = tmp(mxInds(l));
+                        else
+                            mxInds(l) = nan;
+                        end
+                    end
+                    acts{i}(k,:,:) = [mxInds-buffer/2 compInds];
+                    subplot(221);hold all;scatter(mxInds,1:size(tempFull,1),m*5,'g','filled');
+                    multReg = circshift(multReg,[0 -buffer/2]).';
+                    toepReg((k-1)*(buffer+1) + (1:buffer+1),:) = toeplitz(multReg(:),zeros(1,buffer+1)).';
+                end
+                tempFull1 = tempFull.';
+                %wein = (toepReg*toepReg' + lambda*eye((buffer+1)*numel(locs)))\conj(toepReg)*tempFull1(:);
+                [wein,weinInt] = regress(tempFull1(:),toepReg.');%
+                weinInt = diff(weinInt.').'/2;
+                sig = (abs(wein).^2)./(abs(weinInt).^2);
+                wein(sig < 1) = 0;
+                weinOld = weinOld.';
+                oldIm = reshape(weinOld(:).'*toepReg,size(tempFull1));
+                newIm = reshape(wein.'*toepReg,size(tempFull1));
+                [norm(tempFull1(:)-oldIm(:)) norm(tempFull1(:)-newIm(:))]
+            end
+            filts{i} = reshape(wein,[buffer+1 numel(locs)]);
+            subplot(224);sPlot(filts{i}.',[],0);%[t complex(abs(t),angle(t)/6)],[],0);
+            subplot(223);imagesc(complexIm(newIm.'));
+            subplot(222);imagesc(complexIm((tempFull1-newIm).'));drawnow;pause(1);
+        else
+            pause(1);
+        end
+end

hippo/findPlacesNon.m~

@@ -0,0 +1,72 @@
+function [filts cents acts] = findPlacesNon(fields,thresh,sign)
+warning off;
+lambda = 0;
+buffer = 24;%floor(size(fields,3)/2/10);
+if ~exist('sign','var') || sign > 0
+    sign = 1;
+end
+for i = 1:size(fields,1)
+    filts{i} = [];acts{i} = [];cents{i} = [];
+%    for j = 1:2
+        tempFull = sign*squeeze(fields(i,:,:));%(j-1)*size(fields,3)/2+(1:size(fields,3)/2)));
+        tempFull = [zeros(size(tempFull,1),buffer/2) tempFull zeros(size(tempFull,1),buffer/2)];
+        tempFullF = filtfilt(gausswin(buffer/2),sum(gausswin(buffer/2)),tempFull')';
+        tempF = mean(abs(tempFullF));
+        [~,locs] = findpeaks(tempF,'minpeakheight',thresh,'minpeakdistance',buffer/2);
+        subplot(221);hold off;imagesc(complexIm(tempFull));hold all;scatter(locs,ones(numel(locs),1)*size(tempFull,2)/2,'w','filled');
+        if numel(locs)
+            wein = zeros(numel(locs),buffer+1);
+            tempM = mean(tempFull);
+            for k = 1:numel(locs)
+                wein(k,:) = tempM(locs(k)+(-buffer/2:buffer/2));
+            end
+            for m = 1:3
+                if min(size(wein)) == 1
+                    wein = reshape(wein,[ buffer+1 numel(locs)]).';
+                    for k = 1:size(wein,1)
+                        mx = round(sum((1:size(wein,2)).*abs(wein(k,:)))/sum(abs(wein(k,:))));%max(abs(wein(k,:)));
+                        wein(k,:) = circshift(wein(k,:),[0 buffer/2+1-mx]);
+                    end
+                end
+                wein = bsxfun(@rdivide,wein,sqrt(sum(abs(wein).^2,2)));
+                weinOld = wein;
+                toepReg = zeros(numel(locs)*(buffer+1),numel(tempFull));
+                for k = 1:numel(locs)
+                    multReg = zeros(size(tempFull,1),size(tempFull,2));
+                    inds = locs(k)+(-buffer/4:buffer/4);
+                    mxInds = zeros(size(tempFull,1),1);absInds = mxInds;compInds = mxInds;
+                    for l = 1:size(tempFull,1)
+                        tmp = circshift(conv(conj(fliplr(wein(k,:))),tempFull(l,:),'full'),[0 -buffer/2]);
+                        tmp(~ismember(1:numel(tmp),inds)) = 0;
+                        [absInds(l),mxInds(l)] = max(abs(tmp));
+                        if absInds(l) > thresh && mxInds(l) ~= min(inds) && mxInds(l) ~= max(inds)
+                            multReg(l,max(1,mxInds(l))) = tmp(mxInds(l));
+                            compInds(l) = tmp(mxInds(l));
+                        else
+                            mxInds(l) = nan;
+                        end
+                    end
+                    acts{i}(k,:,:) = [mxInds compInds];
+                    subplot(221);hold all;scatter(mxInds,1:size(tempFull,1),m*5,'g','filled');
+                    multReg = circshift(multReg,[0 -buffer/2]).';
+                    toepReg((k-1)*(buffer+1) + (1:buffer+1),:) = toeplitz(multReg(:),zeros(1,buffer+1)).';
+                end
+                tempFull1 = tempFull.';
+                [wein,weinInt] = regress(tempFull1(:),toepReg.');%(toepReg*toepReg' + lambda*eye((buffer+1)*numel(locs)))\conj(toepReg)*tempFull1(:);
+                weinInt = diff(weinInt.').'/2;
+                sig = (abs(wein).^2)./(abs(weinInt).^2);
+                wein(sig < 1) = 0;
+                filts{i} = wein;
+                cents{i}
+                weinOld = weinOld.';
+                oldIm = reshape(weinOld(:).'*toepReg,size(tempFull1));
+                newIm = reshape(wein.'*toepReg,size(tempFull1));
+                [norm(tempFull1(:)-oldIm(:)) norm(tempFull1(:)-newIm(:))]
+            end
+            t = reshape(wein,[buffer+1 numel(locs)]).';
+            subplot(224);sPlot(t,[],0);%[t complex(abs(t),angle(t)/6)],[],0);
+            subplot(223);imagesc(complexIm(newIm.'));
+            subplot(222);imagesc(complexIm((tempFull1-newIm).'));drawnow;
+        end
+%    end
+end

hippo/getData.m

@@ -5,7 +5,7 @@
 % you store the data files.
 
 %data_root = '/media/Expansion Drive/redwood/';%'/media/work/hippocampus/';%
-data_root1 = '/media/Expansion Drive/redwood/';%KenjiMizuseki/';
+data_root1 = '/media/Expansion Drive/KenjiMizuseki/';
 %fname = 'hippo.h5';%'96elec.h5';%
 %padding = 0;%20000;
 %info = hdf5info([data_root fname '.h5']);