cov_matrix.m

Fs=44100;
dr=dir('*.wav');
if length(dr)>0
    nm=[];
    for fc=1:length(dr)
        fln=dr(fc).name;
        lfln=length(fln);
        ns=fln(1:lfln-4);
        nm=[nm str2num(ns)];
    end
    [nms ind]=sort(nm);

    L=length(dr);
    ca=cell(L,1);
    lca=zeros(L,1);
    for fc=1:length(dr)
        [s1 Fs1]=wavread(dr(ind(fc)).name);
        s=resample(s1,2,1);
        ca{fc}=s;
        lca(fc)=length(s);
     
    end
    %soundsc(ca,44100);
    
end

Fs=Fs1*2;

% make equal length specters
%t01=0.1;
%n01=round(t01*Fs);
n01=4096;
n0=8000;
fab=cell(L,1);
for fc=1:L
    fca=fft(ca{fc}(n0:n0+n01-1));
    fabt=abs(fca(1:(n01/2))).^2;
    fab{fc}=fabt/sum(fabt);
    
end



% covariation matrix:
cm=zeros(L,L);

for fc1=1:L
    fab1=fab{fc1};
    for fc2=1:L
        fab2=fab{fc2};
        cm(fc1,fc2)=sum(fab1.*fab2);
    end
end

mx=max(max(cm));
mn=min(min(cm));

cmn=(cm-mn)/(mx-mn);

nd=zeros(L,1); % number of detectd
for fc1=1:L
    %cmn(fc1,:)=cmn(fc1,:)/cmn(fc1,fc1);
    nd(fc1)=length(find(cmn(fc1,:)>=cmn(fc1,fc1)));
end

% ind=find(cmn>0.5);
% cmn(ind)=1;

close('all');
figure;
hi=imagesc(cmn);
axis equal;

%figure;
% plot(nd,'x');
% plot(fab{24},'b-')
% hold on
% plot(fab{31},'r-')

w=ones(n01/2,L); % weights
p=1.3;
r=0.2;
wm=0.2;
for it=1:1000
    for fc=1:L
        ind1=find(p*cmn(fc,:)>=cmn(fc,fc));
        f1=find(ind1~=fc);
        ind=ind1(f1);
        for indc=1:length(ind)
            w(:,fc)=w(:,fc)-r*(fab{ind(indc)}-fab{fc});
            w(:,fc)=(n01/2)*w(:,fc)/sum(w(:,fc));
            w(:,fc)=w(:,fc).*(w(:,fc)>=0);
            
        end

    end
    
    % update matrix:
    for fc1=1:L
        fab1=fab{fc1};
        for fc2=1:L
            fab2=fab{fc2};
            cm(fc1,fc2)=sum(fab1.*w(:,fc1).*fab2);
        end
    end
    mx=max(max(cm));
    mn=min(min(cm));
    cmn=(cm-mn)/(mx-mn);
    
    set(hi,'CData',cmn);
    drawnow;
    
end

cm0=cm;

dg=zeros(L,1);
for fc=1:L
    dg(fc)=cm(fc,fc);
end
mdg=mean(dg);




for fc=1:L
    ka=cm(fc,fc)/mdg;
    w(:,fc)=w(:,fc)/ka;
end

% update matrix:
for fc1=1:L
    fab1=fab{fc1};
    for fc2=1:L
        fab2=fab{fc2};
        cm(fc1,fc2)=sum(fab1.*w(:,fc1).*fab2);
    end
end
mx=max(max(cm));
mn=min(min(cm));
cmn=(cm-mn)/(mx-mn);
imagesc(cmn);


% check
lis=zeros(L,1);
for fc=1:L
        ind1=find((1:L)~=fc);
        mx=max(cmn(fc,ind1));
        lis(fc)=(cmn(fc,fc)-mx)/(cmn(fc,fc)+mx);
end
plot(lis,'x')
mnlis=min(lis)
mxlis=max(lis)

% where w is not 1
wn1=zeros(L,1);
for fc=1:L
    wn1(fc)=std(w(:,fc)-1);
    
end
plot(wn1,'rx');