cov_matrix_minimax_perpendicular.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);
scsm=5;
fab=cell(L,scsm);
fcs=(nst:nen)+49;
scs=1:scsm;
for sc=scs
    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));
        %fabt=real(fca(fmni:fmxi));
        %fabt=fca(fmni:fmxi);
        fab{fc11,sc}=fabt/sqrt(sum(abs(fabt).^2));

        fc11=fc11+1;
    end
end



w=zeros(nfi,L);
for fc=1:L
    fc11=1;
    c=zeros(nfi,scsm*(L-1));
    for sc=scs
        for fc1=1:L
            if fc1~=fc
                c(:,fc11)=fab{fc1,sc};
                fc11=fc11+1;
            end
        end
    end
    M=zeros(scsm*(L-1),scsm*(L-1));
    v=zeros(scsm*(L-1),1);
    %f123=(fab{fc,1}+fab{fc,2}+fab{fc,3})/3;
    f123=fab{fc,1};
    for fc1=1:(scsm*(L-1))
        %v(fc1)=fab{fc,1}'*c(:,fc1);
        v(fc1)=f123'*c(:,fc1);
        for fc2=1:(scsm*(L-1))
            M(fc1,fc2)=c(:,fc1)'*c(:,fc2);
        end
    end
    %r=M\v;
    r=pinv(M)*v;
    sm=zeros(nfi,1);
    for fc1=1:(scsm*(L-1))
        sm=sm+r(fc1)*c(:,fc1);
    end
    %w(:,fc)=fab{fc,1}-sm;
    w(:,fc)=f123-sm;
    w(:,fc)=w(:,fc)/sqrt(w(:,fc)'*w(:,fc));
end

nnn=nst:nen;
fff=linspace(fmn,fmx,size(w,1));

imagesc(fff,nnn,w');
set(gca,'YDir','normal');
xlabel('frequency, Hz');
ylabel('note number');
colorbar;

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