test_w.m

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


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

% normalize w
wn=zeros(size(w));
v1=zeros(1,3);
v1ma=[];
for fc=1:L
    for sc=1:3
        v1(sc)=sum(w(:,fc).*fab{fc,sc});
    end
    v1m=mean(v1);
    wn(:,fc)=w(:,fc)/v1m;
    v1ma=[v1ma v1m];
end



% %check notes with itself
% sp1=3;
% gr=zeros(L,sp1);
% for sp1c=1:sp1
%     for fc=1:L
%         gr(fc,sp1c)=sum(w(:,fc).*fab{fc,sp1c});
%     end
% end
% close('all');
% plot(gr(:,1),'b-');
% hold on;
% plot(gr(:,2),'r-');
% plot(gr(:,3),'g-');
% xlabel('note');
% title('with itself');
% legend('1 piece','2 piece','3 piece');

cm=zeros(L,L,3);
for sc=1:3
    for fc1=1:L
        wt=wn(:,fc1);
        for fc2=1:L
            fab2=fab{fc2,sc};
            cm(fc1,fc2,sc)=sum(wt.*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

% figure;
% imagesc(cmu(:,:,1));
% colorbar;
% 
% figure;
% imagesc(cmu(:,:,2));
% colorbar;
% 
% figure;
% imagesc(cmu(:,:,3));
% colorbar;


dg=[];
for sc=1:3
    for fc=1:L
        dg=[dg cm(fc,fc,sc)];
    end
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 sc=1:3
    for fc1=1:L
        for fc2=1:L
            if fc1~=fc2
                ndg=[ndg cm(fc1,fc2)];
            end
        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) ]);

% close('all');
% plot(dg,'b-');
% hold on;
% plot(ndg,'r-');
save('wn.mat','wn');