cov_matrix_linprog.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

cm0=cm;

n02=n01/2;

% w(1,1)
% w(2,1)
% w(3,1)
% ...
% w(2048,1)
% w(1,2)
% w(2,2)
% w(3,3)


w0=ones(L*n02,1); % initial guess

kd=-1/L; % diag coefficient
knd=2/(L*(L-1)); % not diag coefficient
f=ones(L*n02,1); % initial guess
for fc=1:L
    fp=zeros(n02,1);
    for fc1=1:L
        if fc1==fc
            fp=fp+kd*fab{fc}.*fab{fc1};
        else
            fp=fp+knd*fab{fc}.*fab{fc1};
        end
    end
    f((fc-1)*n02+1:fc*n02)=fp;
end

% kdn*not_diag_element-kd*diag_element -> min


A=zeros(L,L*n02);
for fc=1:L
    A(fc,(fc-1)*n02+1:fc*n02)=fab{fc}.^2;
end
b=ones(L,1);


% x = linprog(f,A,b,Aeq,beq,lb,ub,x0)
lb=zeros(size(w0));
ub=50*ones(size(w0));
options=optimset('Display','iter','MaxIter',100);
[w,fval] = linprog(f,A,b,[],[],lb,ub,w0,options);
%[w,fval] = linprog(f,[],[],[],[],lb,ub,w0,options);


w_lin=w;
w=zeros(n02,L);
for fc=1:L
    w(:,fc)=w_lin((fc-1)*n02+1:fc*n02);
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);
    end
end

imagesc(cm-cm0);

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) ]);