cov_matrix_linprog_sep.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;

uw=100; % up for weights
w=(uw/2)*ones(n02,L);
%options=optimset('Display','iter','MaxIter',1000);
options=optimset('MaxIter',850);
lb=zeros(n02,1);
ub=uw*ones(n02,1);
kd=-1/1; % diag coefficient
knd=1/(L-1); % not diag coefficient
fc0=50;
for fc=1:L
    fc
    f=zeros(n02,1);
    b=0.5*ones(L,1);
    b(fc)=-1;
    A=zeros(L,n02);
    for fc1=1:L
        if fc1==fc
            f=f+kd*fab{fc}.*fab{fc1};
            A(fc,:)=-(fab{fc}.*fab{fc1})';
        else
            f=f+knd*fab{fc}.*fab{fc1};
            A(fc,:)=(fab{fc}.*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);
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);
colorbar;
% close('all');
% plot(cm0(fc0,:),'b-');
% plot(cm(fc0,:),'r-');