cov_matrix_minimax_3pieces.m

global L fc fab ind

nst=-24; % start note
nen=35; % end note

disp(' ');
disp('range:');
disp([number2note(nst) '   ' number2note(nen)]);
disp(['note numbers: ' num2str(nst) '  ' num2str(nen)]);
disp(['note order numbers: ' num2str(nst+49) '  ' num2str(nen+49)]);

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;

% for spectrum cut (not need all):
fmx=440*2^((nen+1)/12); % maximal frequency
fmxi=round(n02*fmx/(Fs/2)); % maximal index
fmn=440*2^((nst-1)/12); % maximal frequency
fmni=round(n02*fmn/(Fs/2)); % minimal index
nfi=fmxi-fmni+1; % number of elements in cuted spectrum

%n0=8000;
n00=100;
L=nen-nst+1;
%fab=cell(L,1);
fab=cell(L,3);
fcs=(nst:nen)+49;
for sc=1:3
    fc11=1;
    for fc=fcs
        i1=n00+1+(sc-1)*n01;
        i2=n00+sc*n01;
        lcat=length(ca{fc});
        if i2>lcat;
            i2=lcat;
            i1=i2-n01+1;
        end
        %catmp=ca{fc}(n0:n0+n01-1);
        catmp=ca{fc}(i1:i2);
        %catmp=ca{fc}(n0:n0+n01-1);
        fca=fft(catmp);
        %fabt=abs(fca(1:(n01/2)));
        fabt=abs(fca(fmni:fmxi));
        fab{fc11,sc}=fabt/sum(fabt);

        fc11=fc11+1;
    end
end



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

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

cm0=cm;

n02=n01/2;

uw=1; % up for weights
%w=(uw/2)*ones(n02,L);
w=(uw/2)*ones(nfi,L);
%options=optimset('Display','iter','MaxIter',1000);
options=optimset('MaxIter',2000);
%options1=optimset('MaxFunEvals',20,'MaxIter',20,'Display','iter');
options1=optimset('MaxIter',40,'Display','iter-detailed');
lb=0*ones(nfi,1);
ub=uw*ones(nfi,1);
% kd=-2*1/1; % diag coefficient
% knd=1/(L-1); % not diag coefficient
% fc0=50;
tic;
for fc=1:L
    
%     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(@mnmxfun1,w(:,fc),[],[],[],[],lb,ub,[],options1);
    jr=fc/L;
    toc1=toc;
    v=jr/toc1; % velosity
    et=(1-jr)/v; % estimated time
    disp(['wait ' num2str(et,'%10.0f') ' s']);
    
end

save('w.mat','w');

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

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

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

r3=ceil(3*rand);
r3n=find((1:3)~=r3);
sc1=r3n(1)
sc2=r3n(2)
cmr=zeros(L,L);
for fc1=1:L
    fab1=fab{fc1,sc1};
    for fc2=1:L
        fab2=fab{fc2,sc2};
        %cm(fc1,fc2)=sum(fab1.*w(:,fc1).*fab2);
        cmr(fc1,fc2)=sum(w(:,fc1).*fab2);
    end
end


cmru=zeros(L,L);
for fc1=1:L
    for fc2=1:L
        if cmr(fc1,fc2)>=cmr(fc1,fc1)
            cmru(fc1,fc2)=1;
        end
    end
end


imagesc(cmru);
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) ]);
%