mfbox_denoise.m

toolbox for spm 5 for data, model free analysis

function [sig]=mfbox_denoise(noisesig,m,k,delayspeed,clustermethod,reducemethod,undelaymethod,sortmethod,choosemethod,gamma,repeat)
%denoise_locca - denoise with local component analysis
%
% Copyright by Peter Gruber and Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the
% GNU GENERAL PUBLIC LICENSE, see gpl.txt
% ------------------------------------------------------------------------------------------------
% REFERENCES:
%
% P. Gruber, K. Stadlthanner, A.M. Tomé, A.R. Teixeira, F.J. Theis, C.G.
% Puntonet, and E.W. Lang. Denoising using local ICA and a generalized
% eigendecomposition with time-delayed signals. In Proc. ICA 2004, volume 3195
% of Lecture Notes in Computer Science, pages 992-1000, Granada, Spain, 2004.
%
% P. Gruber, F.J. Theis, K. Stadlthanner, E.W. Lang, A.M. Tomé, and A.R.
% Teixeira. Denoising using local ICA and kernel-PCA. In Proc. IJCNN 2004, pages
% 2071-2076, Budapest, Hungary, 2004.
%
% P. Gruber, F.J. Theis, A.M. Tomé, and E.W. Lang. Automatic denoising using
% local independent component analysis. In Proc. EIS 2004, Madeira, Portugal,
% 2004.
%

if (nargin<4), delayspeed = 1; end
if (nargin<5), clustermethod = 'kmeans'; end
if (nargin<6), reducemethod = 'pca'; end
if (nargin<7), undelaymethod = 'avg'; end
if (nargin<8), sortmethod = 'var'; end
if (nargin<9), choosemethod = 'mdl'; end
if (nargin<10), gamma = 1; end
if (nargin<11), repeat = 1; end

if (iscell(noisesig))
    grid = noisesig{2};
    scale = ones(size(grid,1),1);
    for i=1:size(grid,1)
        v = unique(grid(i,:));
        scale(i) = 1/min(v(2:end)-v(1:(end-1)));
    end
    grid = diag(scale)*grid;
    noisesig = noisesig{1};
    mi = min(floor(grid),[],2)-1;
    ma = max(ceil(grid),[],2)+1;
    grid = 1+grid-repmat(mi,1,size(grid,2));
    gr = mat2cell(round(grid),ones(1,size(grid,1)),size(grid,2));
    gr = sub2ind(ma'-mi'+1,gr{:});
    newnoisesig = reshape(mfbox_transgrid(grid,noisesig(:), ...
        mfbox_mkgrid(ma'-mi'+1)'),ma'-mi'+1);
    newnoisesig = mfbox_denoise(newnoisesig,m,k,delayspeed,clustermethod, ...
        reducemethod,undelaymethod,sortmethod,choosemethod,gamma,repeat);
    sig = newnoisesig(gr);
    return
end

sigsize = size(noisesig);

orgsigsize = sigsize;

savedsig = delay(noisesig,m,delayspeed);
savedmv = mean(savedsig,2);
for i=1:size(savedsig,1), savedsig(i,:) = savedsig(i,:)-savedmv(i); end

for rcount=1:floor(repeat)
    if (rcount>1)
        rmix = repeat-floor(repeat);
        delayedsig = delay((1-rmix)*noisesig+(rmix)*newsig,m,delayspeed);
        mv = mean(delayedsig,2);
        for i=1:size(delayedsig,1), delayedsig(i,:) = delayedsig(i,:)-mv(i); end
    else
        delayedsig = savedsig; mv = savedmv;
    end

    [delayedsig,comp] = denoise_delayedcoordinates(delayedsig,k, ...
        clustermethod,reducemethod,choosemethod,sortmethod,gamma);

    newsig = undelay(delayedsig+mv*ones(1,size(delayedsig,2)),sigsize, ...
        m,delayspeed,undelaymethod);
end

sig = newsig;

return


function p=choosenumberofcomponents(signal,D,choosemethod,gamma)

[m,samples] = size(signal);

V = sum(triu(D'));

p = 1;
if (m==1), return; end
switch choosemethod
    case 'mdl'
        p = mfbox_mdl(D,samples,gamma);
    case 'err'
        p = m;
        while (sum(V(1,p:m))<gamma*sum(V(1,:))) & (p>0)
            p = p-1;
        end;
    case 'abs'
        p = min(ceil(abs(gamma)),m);
    case 'num'
        p = min(max(abs(gamma),0),m);
    case 'cluster'
        if (gamma>0), [ce,cl] = kmeans_cluster(V(1,:),2); p = sum(cl==cl(1));
        else [ce,cl] = cluster(V,2); p = sum(cl==cl(1)); end
end

return


function delaysig=delay(sig,lag,step)

d = size(sig);

if ((length(d)==2) && ((d(1)==1)||(d(2)==1)))
    delaysig = [];
    for i=-lag:lag
        delaysig = [delaysig;circshift(sig',-round(i*step))'];
    end
    return
end

ml = ceil(lag)+2;
linoffs = mfbox_mkgrid(1+(ones(1,length(d))*2*ml))'-ml;
linoffs = linoffs(:,sum(linoffs.^2)<=(lag^2));

lo = size(linoffs,2);
%linoffs = ceil(step*linoffs);
linoffs = ceil(step*(linoffs.*repmat(sum(linoffs.^2)/lag^2,length(d),1)));

if (lo==0)
    delaysig = []; return
end

nd = prod(d);

delaysig = zeros(nd,lo);
t = circshift(sig,-linoffs(:,1));
delaysig(:,1) = t(:);
for j=2:lo
    t = circshift(t,linoffs(:,[j-1:j])*[1;-1]);
    delaysig(:,j) = t(:);
end

delaysig = delaysig';

return


function sig=delocalize(localsig,permutation)

[dummy,p] = sort(permutation);
sig = localsig(:,p);

return


function [sig,rk,components]=denoise_delayedcoordinates(noisesig,k,clustermethod,reducemethod,choosemethod,sortmethod,gamma)

[m,samples] = size(noisesig);

meanvalues = mean(noisesig,2);
for i=1:size(noisesig,1), noisesig(i,:) = noisesig(i,:)-meanvalues(i); end

[localsig,parts,permut] = localize(noisesig,k,m,clustermethod);

ncomp = size(parts,1);

rk = zeros(1,ncomp);

for i=1:ncomp
    lsig = localsig(:,parts(i,1):parts(i,2));
    lmean = mean(lsig,2);
    for j=1:size(lsig,1), lsig(j,:) = lsig(j,:)-lmean(j); end

    [V,W,D] = sortcomponents(lsig,reducemethod,sortmethod);
    p = choosenumberofcomponents(V*lsig,D,choosemethod,gamma);

    r = size(V,1);
    R = eye(r);
    R(:,p+1:r) = 0;

    rk(i) = p;
    lsig = W*R*V*lsig;

    localsig(:,parts(i,1):parts(i,2)) = lsig+lmean*ones(1,size(lsig,2));
end

sig = delocalize(localsig,permut)+meanvalues*ones(1,samples);

return


function [localdata,parts,permutation]=localize(sig,k,m,method)

[dim,samples] = size(sig);

if (k==1), localdata = sig; parts = [1,samples]; permutation = 1:samples; return; end

if (m>samples), m = samples; end

switch method
    case 'kmeans'
        [cent,cl] = cluster_kmeans(sig,k);
    case 'fast'
        numkmeans = max(max(500,ceil(samples/500)),n);
        rsig = sig(:,randperm(numkmeans));
        [cent,cl] = cluster_kmeans(rsig,n);
        [err,cl] = min(repmat(sum(cent.^2,1)',1,samples)-2*cent'*sig);
    case 'seq'
        [cent,cl] = cluster(sig,k);
end

p = zeros(size(cent,2),1);
for i=1:size(cent,2) p(i) = sum(cl==i); end

while ((min(p)<m) || (sum(isfinite(p))>k))
    [i,im] = min(p);
    q = zeros(size(cent,2),1);
    for i=1:size(cent,2)
        if (isfinite(p(i))), q(i) = dist(cent(:,im),cent(:,i),2);
        else q(i) = +Inf;
        end
    end
    [i,j] = sort(q);
    cl = cl+(cl==j(1))*(j(2)-j(1));
    p(j(2)) = p(j(2))+p(j(1));
    p(j(1)) = +Inf;
end

k = sum(isfinite(p));
pt = zeros(k,3);
[i,j] = sort(p);
pt(:,3) = cumsum(i(1:k));
pt(:,1) = j(1:k);
pt(1,2) = 1; pt(2:k,2) = pt(1:k-1,3)+1;

i = zeros(samples,1);
for j=1:k, i = i+(cl==pt(j,1))*(j-pt(j,1)); end
cl = cl+i; clear i

localdata = zeros(dim,samples);
permutation = zeros(samples,1);

nn = 1:samples;
for i=1:k
    localdata(:,pt(i,2):pt(i,3)) = sig(:,cl==i);
    permutation(pt(i,2):pt(i,3)) = nn(cl==i);
end

parts = pt(:,2:3);

return


function [V,W,D]=sortcomponents(signal,reducemethod,ratingmethod)

[m,samples] = size(signal);

switch reducemethod
    case 'pca'
        [W,D] = eig(cov(signal'));
        V = W';
        ksignal = V*signal;
    case 'mpencil'
        [W,D] = eig(cov(signal'));
        E = D; E(E>0) = 1./sqrt(E(E>0));
        V = W*E*W';
        c = max(10,ceil(samples/1000));
        c = exp(-[-c:1:c].^2)/sqrt(2*pi);
        [P,Q] = eig(cov(matconv(V*signal,c)'));
        E = D; E(E>0) = sqrt(E(E>0));
        W = W*E*W'*P;
        W = W./(ones(m,1)*sqrt(sum(W.^2)));
        W = ort(W);
        V = W^-1;
        ksignal = V*signal;
    case 'pearsonica'
        [ksignal,W,V] = mfbox_fasticapearson(signal',0.0001,200,[2.6,4],[0,1], ...
            size(signal,1),1,size(signal,1),'off');
        W = W./(ones(m,1)*sqrt(sum(W.^2)));
        W = ort(W);
        V = W^-1;
        ksignal = V*signal;
    case 'mdlica'
        [E,D] = eig(cov(signal'));
        D = sum(D,1); D = D(m:-1:1); E = E(:,m:-1:1);
        p = mfbox_mdl(D,samples,1);
        ksignal = E'*signal;
        if (p==1)
            ksignal = ksignal(1,:);
            V = 1; W = 1;
        else
            [ksignal,W,V] = mfbox_fasticapearson(ksignal(1:p,:)',0.0001, ...
                200,[2.6,4],[0,1],p,1,p,'off');
        end
        V = [V,zeros(p,m-p)]*E';
        W = E*[W;zeros(p,m-p)];
        W = W./(ones(m,1)*sqrt(sum(W.^2)));
        V = W^-1;
        ksignal = V*signal;
end

k = size(ksignal,1);
R = eye(k);
R = R(k:-1:1,:);

switch ratingmethod
    case 'kurt'
        [dummy,kp] = sort(kurtosis(ksignal').^2);
        [dummy2,kp] = sort(kp);
    case 'maxkurt'
        [dummy,kp] = sort(-kurtosis(ksignal'));
        [dummy2,kp] = sort(kp);
    case 'minkurt'
        [dummy,kp] = sort(kurtosis(ksignal'));
        [dummy2,kp] = sort(kp);
    case 'comb'
        [dummy,kp] = sort((var(ksignal').*kurtosis(ksignal')).^2);
        [dummy2,kp] = sort(kp);
    case 'auto'
        [dummy,kp] = sort(var(autocor(ksignal')));
        [dummy2,kp] = sort(kp);
    otherwise
        [dummy,kp] = sort(var(ksignal'));
        [dummy2,kp] = sort(kp);
end
R = R(:,kp);
V = R*V;
W = W*R';

switch reducemethod
    case 'pca'
        D = ones(m,1)*var((R*ksignal)');
    otherwise
        T = W;
        C = cov((R*ksignal)');
        D(1,:) = diag(C)';
        for i=1:m
            S = ort(V(i+1:m,:));
            S = ort(S(m-i+1:m,:));
            d = S*T*C*T'*S';
            D(i,1:i) = eig(d(1:i,1:i))';
            D(i,i+1:m) = eig(d(i+1:m,i+1:m))';
        end
end

return


function N=matconv(M,C)

[d,m] = size(M);
c = length(C);
N = zeros(d,m+c-1);
for i=1:size(M,1)
    N(i,:) = conv(M(i,:),C(:));
end

return


function sig=undelay(delaysig,d,lag,step,recon)

if (length(d)==2)
    if ((d(1)==1)||(d(2)==1))
        q = lag^2+1-(-lag:lag).^2; q = (2*lag+1)*q/sum(q);
        for i=-lag:lag
            j = i+lag+1;
            b = round(i*step);
            switch recon
                case 'avg'
                    delaysig(j,:) = circshift(delaysig(j,:)',b)';
                case 'quad'
                    delaysig(j,:) = q(j)*circshift(delaysig(j,:)',b)';
            end
        end
        sig = sum(delaysig,1)/size(delaysig,1);
        return
    end
end

pot = 0;
switch recon
    case 'avg'
        pot = 0;
    case 'quad'
        pot = 2;
    otherwise
        pot = 0;
end

ml = ceil(lag)+2;
linoffs = mfbox_mkgrid(1+(ones(1,length(d))*2*ml))'-ml;
weights = sum(linoffs.^2);
linoffs = linoffs(:,weights<=(lag^2));
weights = sqrt(weights(weights<=(lag^2))).^pot;

lo = size(linoffs,2);
%linoffs = ceil(step*linoffs);
linoffs = ceil(step*(linoffs.*repmat(sum(linoffs.^2)/lag^2,length(d),1)));

if (lo==0)
    delaysig = []; return
end

nd = prod(d);
ld = length(d);

weights = 1+weights-min(weights);
weights = weights/sum(weights);

sig = zeros(d);
delaysig = delaysig';

for k=1:lo
    t = reshape(delaysig(:,k),d);
    t = circshift(t,linoffs(:,k));
    sig = sig+(weights(k)*t);
end

return


function [centroids,clusters,err]=cluster(D,n)

epochs = 100;

[dim,samples] = size(D);

if (n==1)
    centroids = sum(D,2)/samples;
    clusters = ones(samples,1);
    err = sum(sqrt(sum((D-centroids*ones(1,samples)).^2,1)))/samples;
    return
end

partition = zeros(1,n+1); centers = zeros(1,n);

for i=2:n, partition(i) = floor(samples/n)*(i-1); end
partition(n+1) = samples;

for t=1:epochs
    oldpart = partition;
    for i=1:n, centers(i) = partition(i)+ceil((partition(i+1)-partition(i))/2); end
    for i=1:n-1
        num = centers(i+1)-centers(i)+1;
        if (num==2), partition(i+1) = centers(i);
        else
            dists = zeros(1,num);
            for j=1:num-1
                dists(j+1) = dists(j)+ ...
                    sqrt(sum((D(:,centers(i)+j-1)-D(:,centers(i)+j)).^2));
            end
            partition(i+1) = centers(i)+sum(dists<(dists(num)/2));
        end
    end
    if all(partition==oldpart) break; end
end

centroids = zeros(dim,n); clusters = zeros(samples,1); err = 0;
for i=2:n+1
    s = partition(i-1)+1;
    e = partition(i); n = e-s+1;
    centroids(:,i-1) = sum(D(:,s:e),2)/n;
    clusters(s:e) = ones(1,n)*(i-1);
    err = err+sum(sqrt(sum((D(:,s:e)-centroids(:,i-1)*ones(1,n)).^2,1)));
end

err = err/samples;

return


function [centroids,clusters]=cluster_kmeans(data,n)

[dim dlen] = size(data);
n = min(dlen,n);

temp = randperm(dlen);
centroids = data(:,temp(1:n));

clusters = zeros(1,dlen);

iter = 0;
old_clusters = zeros(1,dlen);
while (iter<100)
    [dummy clusters] = min(repmat(sum((centroids.^2),1)',1,dlen)-2*centroids'*data);
    for i=1:n,
        f = find(clusters==i);
        s = length(f);
        if (s), centroids(:,i) = sum(data(:,f),2)/s; end
    end
    if ((iter>0)&&all(old_clusters==clusters))
        break;
    end
    old_clusters = clusters;
    iter = iter+1;
end

[qerrs clusters] = min(repmat(sum((centroids.^2),1)',1,dlen)-2*centroids'*data);
clusters = clusters(:);

return


function A=ort(A)

[rows,columns] = size(A);

if (rows>columns)
    A = ort(A')';
    return
end

for i=1:rows
    A(i,:) = A(i,:)/norm(A(i,:));
end

while (rows<columns)
    m = zeros(rows+1);
    m(1:rows,1:rows+1) = A(1:rows,1:rows+1);
    n = zeros(1,columns);
    for i=1:rows+1
        m(rows+1,i) = 1;
        n(1,i) = det(m);
        m(rows+1,i) = 0;
    end
    if (norm(n)>0), n = n/norm(n); end
    A = [A;n];
    rows = rows+1;
end

return


function err=dist(A,B,pot)

A = reshape(A,1,prod(size(A)));
B = reshape(B,1,prod(size(B)));

err = sum(sum(sqrt((A-B).*conj(A-B)).^pot).^(1/pot))/size(A,2);

return


function retval=kurtosis(x)

[nr,nc] = size(x);
x = x-ones(nr,1)*mean(x);
retval = zeros(1,nc);
s = std(x);
ind = find(s>0);
retval(ind) = sum(x(:,ind).^4)./(nr*s(ind).^4)-3;

return