cov_matrix_3.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)));
    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



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(cm);
axis equal;

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

w=ones(n01/2,L); % weights
p=1.5;
r=0.3;
r1=0.3;
wm=0.2;
for it=1:100
    for fc=1:L
        ind1=find(p*cm(fc,:)>=cm(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
        
    set(hi,'CData',cm);
    drawnow;
    
end
cm0=cm;
w0=w;

cm=cm0;
w=w0;

dg=zeros(L,1);
for fc=1:L
    dg(fc)=cm(fc,fc);
end
%plot(dg)
for fc=1:L
    w(:,fc)=w(:,fc)/dg(fc);
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
    
    dg=zeros(L,1);
    
for fc=1:L
    dg(fc)=cm(fc,fc);
end
%plot(dg)
set(hi,'CData',cm);
drawnow;

max(max(cm))