cov_matrix_minimax.m

global L fc fab ind

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;
n02=n01/2;
fmx=4000; % maximal frequency
fmxi=round(n02*fmx/(Fs/2)); % maximal
n0=8000;
fab=cell(L,1);
for fc=1:L
    fca=fft(ca{fc}(n0:n0+n01-1));
    %fabt=abs(fca(1:(n01/2)));
    fabt=abs(fca(1:fmxi));
    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

cm0=cm;

n02=n01/2;

uw=1; % up for weights
%w=(uw/2)*ones(n02,L);
w=(uw/2)*ones(fmxi,L);
%options=optimset('Display','iter','MaxIter',1000);
options=optimset('MaxIter',2000);
%options1=optimset('MaxFunEvals',20,'MaxIter',20,'Display','iter');
options1=optimset('MaxIter',10,'Display','iter-detailed');
lb=0*ones(fmxi,1);
ub=uw*ones(fmxi,1);
kd=-2*1/1; % diag coefficient
knd=1/(L-1); % not diag coefficient
fc0=50;
for fc=1:L
    fc
    f=zeros(fmxi,1);
    for fc1=1:L
        if fc1==fc
            %f=f+kd*fab{fc}.*fab{fc1};
            f=f+kd*fab{fc1};

        else
            %f=f+knd*fab{fc}.*fab{fc1};
            f=f+knd*fab{fc1};

        end
    end
    
    %[w(:,fc),fval] = linprog(f,A,b,[],[],lb,ub,w(:,fc),options);
    %[w(:,fc),fval] = linprog(f,A,b);
    [w(:,fc),fval] = linprog(f,[],[],[],[],lb,ub,w(:,fc),options);
    %x = fminimax(fun,x0)
    % x = fminimax(fun,x,A,b,Aeq,beq,lb,ub)
    ind=find((1:L)~=fc);
    w(:,fc)=fminimax(@mnmxfun,w(:,fc),[],[],[],[],lb,ub,[],options1);
    
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.*w(:,fc1).*fab2);
        cm(fc1,fc2)=sum(w(:,fc1).*fab2);
    end
end

cmu=zeros(L,L);
for fc1=1:L
    for fc2=1:L
        if cm(fc1,fc2)>=cm(fc1,fc1)
            cmu(fc1,fc2)=1;
        end
    end
end
        

imagesc(cm);
colorbar;
% close('all');
% plot(cm0(fc0,:),'b-');
% plot(cm(fc0,:),'r-');

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

dgm=mean(dg);
dgs=std(dg);
dgmx=max(dg);
dgmn=min(dg);

disp(' ');
disp(' ');
disp(['mean of diagoanal: ' num2str(dgm) ]);
disp(['std of diagoanal: ' num2str(dgs) ]);
disp(['max of diagoanal: ' num2str(dgmx) ]);
disp(['min of diagoanal: ' num2str(dgmn) ]);

ndg=[];
for fc1=1:L
    for fc2=1:L
        if fc1~=fc2
            ndg=[ndg cm(fc1,fc2)];
        end
    end
end

ndgm=mean(ndg);
ndgs=std(ndg);
ndgmx=max(ndg);
ndgmn=min(ndg);

disp(' ');
disp(['mean of not diagoanal: ' num2str(ndgm) ]);
disp(['std of not diagoanal: ' num2str(ndgs) ]);
disp(['max of not diagoanal: ' num2str(ndgmx) ]);
disp(['min of not diagoanal: ' num2str(ndgmn) ]);