cov_matrix_minimax_perpendicular_improve.m

global L fc fab ind

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

load('sma.mat');

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



w=zeros(nfi,L);
for fc=1:L
    fc11=1;
    c=zeros(nfi,3*L-3);
    for sc=1:3
        for fc1=1:L
            if fc1~=fc
                c(:,fc11)=fab{fc1,sc};
                fc11=fc11+1;
            end
        end
    end
    
    nta=false; % if need to add
    for smc=1:length(sma)-2 % 2 simul, 3 simul, 4 simul
        smt=sma{smc};
        fn=find(smt(:,1)==fc);
        lfn=length(fn);
        if lfn>0 % if any notes
            nta=true;
            for fnc=1:lfn % for all notes that satisfy
                st=zeros(n01,1);
                for snt=2:(smc+1)
                    %nst+49 -start note
                    %n0=n0min+(n0max-n0min)*rand;
                    %n0=round(n0);
                    st=st+ca{smt(fn(fnc),snt)+nst+49-1}(n0:n0+n01-1);
                end
                
                fca=fft(st);
                fabt=abs(fca(fmni:fmxi));
                fabt=fabt/sum(fabt);
                c=[c fabt];
                
            end              
            
            
        end
        
    end
    
    if nta % if need to extend
        lc=length(c(1,:));
        M=zeros(lc,lc);
        v=zeros(lc,1);
        %f123=(fab{fc,1}+fab{fc,2}+fab{fc,3})/3;
        f123=fab{fc,1};
        for fc1=1:lc
            %v(fc1)=fab{fc,1}'*c(:,fc1);
            v(fc1)=f123'*c(:,fc1);
            for fc2=1:lc
                M(fc1,fc2)=c(:,fc1)'*c(:,fc2);
            end
        end
        %r=M\v;
        r=pinv(M)*v;
        sm=zeros(nfi,1);
        for fc1=1:lc
            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));
        
    else
        
        M=zeros(3*L-3,3*L-3);
        v=zeros(3*L-3,1);
        %f123=(fab{fc,1}+fab{fc,2}+fab{fc,3})/3;
        f123=fab{fc,1};
        for fc1=1:(3*L-3)
            %v(fc1)=fab{fc,1}'*c(:,fc1);
            v(fc1)=f123'*c(:,fc1);
            for fc2=1:(3*L-3)
                M(fc1,fc2)=c(:,fc1)'*c(:,fc2);
            end
        end
        %r=M\v;
        r=pinv(M)*v;
        sm=zeros(nfi,1);
        for fc1=1:(3*L-3)
            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
end

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