S-MSVQ_Codec2.m

%---------------------------------------------------------------------------
% This script:
% * loads wave file with the speech corpus
% * extracts LPC coefficients for every 20 ms frame (with 10 ms overlap)
% * converts LPCs to LSFs and builds a training codebook: Nx10 matrix of floats
% * calculates coarse codebook plus additional fine codebooks for split MSVQ
% * dumps codebook data to a text file (codebooks.txt)
%
% In case of convergence error - re-run the script.
%
% Wojciech Kaczmarski, SP5WWP
% M17 Project (m17project.org)
% 22/01/2022
%---------------------------------------------------------------------------

%samples per 20 ms frame
%8000*0.02
SPF=160;

%codebook sizes
Q0_size=64;
Q1_size=16;
Q2_size=16;
Q3_size=16;

%draw fancy plots?
b_plot=1;

%Fs has to be 8 kHz
[sample, Fs]=audioread("corpus.wav");
if Fs~=8000
    fprintf("Sample rate has to be 8000 Hz.\nAborting.\n");
    return;
end
len=length(sample);

%window
w=hamming(SPF);

frame=zeros(SPF);
a=zeros(11);
lsf=zeros(10, (len-mod(len, SPF)-SPF)/(SPF/2));

%iterate through the whole sample, but omit last few samples
%samples overlap! thats why the step is set to SPF/2
for i=1:SPF/2:len-mod(len, SPF)-SPF
    frame=sample(i:i+SPF-1).*w;
    [a, g]=lpc(frame, 10);
    lsf(:,round(i/(SPF/2)+1)) = cos(poly2lsf(a)); %cosine domain
end

lsf=lsf.';
siz=size(lsf,1); %training set's size

%vector lists
q1=lsf(:,1:3);
q2=lsf(:,4:6);
q3=lsf(:,7:10);

%------------------Q0------------------
%generate "coarse" vector codebook Q0 based on the input data
[~, Q0]=kmeans(lsf, Q0_size);

%---------------first part----------------
%------------------Qf1_1------------------
sQf1_1=zeros(siz,3);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Q0(:,1:3), lsf(i,1:3));
    sQf1_1(i,1:3)=lsf(i,1:3)-Q0(k,1:3);
end

[~, Qf1_1]=kmeans(sQf1_1, Q1_size);

%------------------Qf1_2------------------
sQf1_2=zeros(siz,3);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Qf1_1(:,1:3), sQf1_1(i,1:3));
    sQf1_2(i,1:3)=sQf1_1(i,1:3)-Qf1_1(k,1:3);
end

[~, Qf1_2]=kmeans(sQf1_2, Q1_size);

%------------------Qf1_3------------------
sQf1_3=zeros(siz,3);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Qf1_2(:,1:3), sQf1_2(i,1:3));
    sQf1_3(i,1:3)=sQf1_2(i,1:3)-Qf1_2(k,1:3);
end

[~, Qf1_3]=kmeans(sQf1_3, Q1_size);

%test - clusters should diverge to origin
if b_plot
    figure(1);
    plot(sQf1_1(:,1), sQf1_1(:,2), '. red');
    hold on;
    plot(sQf1_2(:,1), sQf1_2(:,2), '. green');
    hold on;
    plot(sQf1_3(:,1), sQf1_3(:,2), '. blue');
end

%---------------second part---------------
%------------------Qf2_1------------------
sQf2_1=zeros(siz,3);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Q0(:,4:6), lsf(i,4:6));
    sQf2_1(i,1:3)=lsf(i,4:6)-Q0(k,4:6);
end

[~, Qf2_1]=kmeans(sQf2_1, Q2_size);

%------------------Qf2_2------------------
sQf2_2=zeros(siz,3);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Qf2_1(:,1:3), sQf2_1(i,1:3));
    sQf2_2(i,1:3)=sQf2_1(i,1:3)-Qf2_1(k,1:3);
end

[~, Qf2_2]=kmeans(sQf2_2, Q2_size);

%------------------Qf2_3------------------
sQf2_3=zeros(siz,3);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Qf2_2(:,1:3), sQf2_2(i,1:3));
    sQf2_3(i,1:3)=sQf2_2(i,1:3)-Qf2_2(k,1:3);
end

[~, Qf2_3]=kmeans(sQf2_3, Q2_size);

%------------------Qf2_4------------------
sQf2_4=zeros(siz,3);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Qf2_3(:,1:3), sQf2_3(i,1:3));
    sQf2_4(i,1:3)=sQf2_3(i,1:3)-Qf2_3(k,1:3);
end

[~, Qf2_4]=kmeans(sQf2_4, Q2_size);

%test - clusters should converge to origin
if b_plot
    figure(2);
    plot(sQf2_1(:,1), sQf2_1(:,2), '. red');
    hold on;
    plot(sQf2_2(:,1), sQf2_2(:,2), '. green');
    hold on;
    plot(sQf2_3(:,1), sQf2_3(:,2), '. blue');
    hold on;
    plot(sQf2_4(:,1), sQf2_4(:,2), '. cyan');
end

%---------------third part----------------
%------------------Qf1_1------------------
sQf3_1=zeros(siz,4);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Q0(:,7:10), lsf(i,7:10));
    sQf3_1(i,1:4)=lsf(i,7:10)-Q0(k,7:10);
end

[~, Qf3_1]=kmeans(sQf3_1, Q3_size);

%------------------Qf1_2------------------
sQf3_2=zeros(siz,4);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Qf3_1(:,1:4), sQf3_1(i,1:4));
    sQf3_2(i,1:4)=sQf3_1(i,1:4)-Qf3_1(k,1:4);
end

[~, Qf3_2]=kmeans(sQf3_2, Q3_size);

%------------------Qf1_3------------------
sQf3_3=zeros(siz,4);
for i=1:siz
    %find the closest codebook entry for each vector
    [k,dist] = dsearchn(Qf3_2(:,1:4), sQf3_2(i,1:4));
    sQf3_3(i,1:4)=sQf3_2(i,1:4)-Qf3_2(k,1:4);
end

[idx, Qf3_3]=kmeans(sQf3_3, Q3_size);

%test - clusters should diverge to origin
if b_plot
    figure(3);
    plot(sQf3_1(:,1), sQf3_1(:,2), '. red');
    hold on;
    plot(sQf3_2(:,1), sQf3_2(:,2), '. green');
    hold on;
    plot(sQf3_3(:,1), sQf3_3(:,2), '. blue');
end

%dump it all to a .txt file
format long;
fp=fopen('codebooks.txt', 'w');

%-------------------Q0-------------------
fprintf(fp, "Q0:\n");
for i=1:Q0_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff, %1.12ff, %1.12ff, %1.12ff, %1.12ff, %1.12ff, %1.12ff, %1.12ff},\n", ...
        Q0(i,1), Q0(i,2), Q0(i,3), ...
        Q0(i,4), Q0(i,5), Q0(i,6), ...
        Q0(i,7), Q0(i,8), Q0(i,9), ...
        Q0(i,10));
end
fprintf(fp, "\n");

%------------------Qf1_1-----------------
fprintf(fp, "Qf1_1:\n");
for i=1:Q1_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf1_1(i,1), Qf1_1(i,2), Qf1_1(i,3));
end
fprintf(fp, "\n");

%------------------Qf1_2-----------------
fprintf(fp, "Qf1_2:\n");
for i=1:Q1_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf1_1(i,2), Qf1_2(i,2), Qf1_2(i,3));
end
fprintf(fp, "\n");

%------------------Qf1_3-----------------
fprintf(fp, "Qf1_3:\n");
for i=1:Q1_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf1_1(i,3), Qf1_3(i,2), Qf1_3(i,3));
end
fprintf(fp, "\n");

%------------------Qf2_1-----------------
fprintf(fp, "Qf2_1:\n");
for i=1:Q2_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf2_1(i,1), Qf2_1(i,2), Qf2_1(i,3));
end
fprintf(fp, "\n");

%------------------Qf2_2-----------------
fprintf(fp, "Qf2_2:\n");
for i=1:Q2_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf2_2(i,1), Qf2_2(i,2), Qf2_2(i,3));
end
fprintf(fp, "\n");

%------------------Qf2_3-----------------
fprintf(fp, "Qf2_3:\n");
for i=1:Q2_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf2_3(i,1), Qf2_3(i,2), Qf2_3(i,3));
end
fprintf(fp, "\n");

%------------------Qf2_4-----------------
fprintf(fp, "Qf2_4:\n");
for i=1:Q2_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf2_4(i,1), Qf2_4(i,2), Qf2_4(i,3));
end
fprintf(fp, "\n");

%------------------Qf3_1-----------------
fprintf(fp, "Qf3_1:\n");
for i=1:Q3_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf3_1(i,1), Qf3_1(i,2), Qf3_1(i,3), Qf3_1(i,4));
end
fprintf(fp, "\n");

%------------------Qf3_2-----------------
fprintf(fp, "Qf3_2:\n");
for i=1:Q3_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf3_2(i,1), Qf3_2(i,2), Qf3_2(i,3), Qf3_2(i,4));
end
fprintf(fp, "\n");

%------------------Qf3_3-----------------
fprintf(fp, "Qf3_3:\n");
for i=1:Q3_size
    fprintf(fp, "\t{%1.12ff, %1.12ff, %1.12ff, %1.12ff},\n", ...
        Qf3_3(i,1), Qf3_3(i,2), Qf3_3(i,3), Qf3_3(i,4));
end
fprintf(fp, "\n");

fclose(fp);