clear all

%Volterra speech filter using covariance extension of linear LMS filter for
%audio compression utilizing periodic delay for long term prediction

%overall goal of the design is to measure the energy of the speech signal in blocks of .02*Fs and for various delay values.
%the energy of these delay blocks are compared in order to find which delay produces
%the minimum energy. This delay block is selected and used to generate a predicted signal
%if the predicted signal has a lower energy of the original signal and the signal generated
%by a linear filter in an acceptable amount of time, then adequate signal
%compression is achieved.

[input_signal,Fs] = wavread('female_german');   %load audio file @8kHz

   %Define filter parameters
   
window = Fs*.02;                                %window size is determined by the number of samples taken in 20 milliseconds
reps = floor(length(input_signal)/window);      %determine number of repitions required to evaluate the entire signal

s = input_signal(1:reps*window);                   %signal with adjusted length
M = 3;                                          %order of the filter
D = 150;                                        %Maximum delay of the signal

   %energy parameters
Eplot_s = s(161:length(s)).*s(161:length(s));
Energy_s = sum(Eplot_s);
error = [];
tdelay = [];
Eplot = [];

tic
for c = 1:reps-1;                               %repeat script to cover the entire length of the signal in blocks of length(window) samples

   %construct R covariance matrices
e = [];
sample = s(1+window*(c-1):window*(c+1));        %sample portion of the file
N = length(sample);                             %number of samples being analyzed

   %the R matrices are derived using the Covariance method of the Volterra
   %series expansion formulas
for d = 10:D
R2 = zeros(M+1,M+1);
R3 = zeros(M+1,M+1,M+1);
R4 = zeros(M+1,M+1,M+1,M+1);
R3b = zeros(1,M+1,M+1);

for i = 1+d:M+1+d
   for l = 1+d:M+1+d
       R2(i-d,l-d) = sum(sample(1+(M+1)+N/2-i:(N)-i).*sample(1+(M+1)+N/2-l:(N)-l));
   end
end

for k = 1+d:M+1+d
   for j = k:M+1+d
       for l = 1+d:M+1+d
           R3(k-d,j-d,l-d) = sum(sample(1+(M+1)+N/2-k:(N)-k).*sample(1+(M+1)+N/2-j:(N)-j).*sample(1+(M+1)+N/2-l:(N)-l));
       end
   end
end

for k = 1+d:M+1+d
   for j = k:M+1+d
       for p = 1+d:M+1+d
           for q = p:M+1+d
               R4(k-d,j-d,p-d,q-d) = sum(sample(1+(M+1)+N/2-k:(N)-k).*sample(1+(M+1)+N/2-j:(N)-j).*sample(1+(M+1)+N/2-p:(N)-p).*sample(1+(M+1)+N/2-q:(N)-q));
           end
       end
   end
end        


   %Construct Q quadrants for R super matrix. Note: when constructing matrices
   %the first value is ommited. eg R2(i,j) --> R(0,:) and R(:,0) are not used
   %for constructing the super matrix. furthermore Q2 = Q3'

Q1 = R2(2:M+1,2:M+1);
Q2 = [];

Q4 = [];
Q4vert = [];


for x = 2:M+1
   q2 = [];    
   
   switch x
   
       case M+1                                %m = sum(i) for i = 1:M-1
   q2(:,:) = R3(x,x,2:M+1);                    %switch statements required to retain column/row form for concatenation
   Q2 = horzcat(Q2,q2);            
   
       otherwise
   q2(:,:) = R3(x,x:M+1,2:M+1);    
   Q2 = horzcat(Q2,q2');                       %concatenate q2 to form Quandrant 2 matrix. size is M x m
   end
   
   for y = 2:M+1
       
       if (y == M+1)|| (x == M+1)              %if else statements in place to retain orientation for concatenation
           q4 = [];
       q4(:,:)=R4(x,x:M+1,y,y:M+1);
       Q4vert = vertcat(Q4vert,q4);
       
       else
       q4 = [];
       q4(:,:)=R4(x,x:M+1,y,y:M+1);    
       Q4vert = vertcat(Q4vert,q4');  
       end
   end
   
   
   Q4 = horzcat(Q4,Q4vert);                    %horizontally concatenate Q4horz M times. this creates the Quadrant 4 matrix of size (M^2) x (M^2)
   Q4vert = [];
end

Rtop = horzcat(Q1,Q2);                          %concatenate first two quadrants, creating the top half of super matrix R
Rbot = horzcat(Q2',Q4);                         %concatenate last two quadrants, creating the bottom half of super matrix R
R = triu(vertcat(Rtop,Rbot));                   %concatenate both halves to form Super matrix R

   %Construct b vector


for i = 1+d:M+1+d
   R2b(i-d,1) = sum(sample(1+(M+1)+N/2-i:(N+1)-i).*sample((M+1)+N/2:(N)));
end



for j = 1+d:M+1+d
   for l = 1+d:M+1+d
       R3b(1,j-d,l-d) = sum(sample((M+1)+N/2:(N)).*sample(1+(M+1)+N/2-j:(N+1)-j).*sample(1+(M+1)+N/2-l:(N+1)-l));
   end
end
   
   
b2 = [];
b1 = R2b(2:M+1,1);                              %first portion of b vector is equal to R2(0,:) and is a column vector

for i = 2:M+1
   for j = i:M+1
   b2_R3(:,:) = R3b(1,i,j);                    %Define subvector such that b2 will retain Mx1 size
   b2 = vertcat(b2,b2_R3);                     %concatenate to construct the rest of the b2 subvector which is equal to R3(0,1,:) ... R3(0,M,:)
   end
end

b = vertcat(b1,b2);                             %concatenate R20 and R30 portion of b vector

   %system is in Ax=b format where A = R, x = h, and b = R0 vectors

h = R\b;                                        %construct h vector


h1 = h(1:M);
hid = M+1;
h2 = zeros(M,M);

for i = 1:M
   for j = i:M
       h2(i,j) = h(hid);                       %split h vector into linear and quadratic sections to derive predicted sample
       hid = hid+1;
   end
end


ps2 = zeros(N+M,1);
ps1 = zeros(N/2,1);

   %perform convolution of h with s to obtain the predicted sample

csample = padarray(sample(N/2+1:N),M,'pre');    %pad convolution sample to account for first M samples

for n = M+1:N/2+M
   ps1(n,1) = sum(h1(1:M).*csample(n-1:-1:n-M)); %linear predicted sample
end

for n = M+1:N/2+M
   for i = 2:M+1
       for j = i:M+1
           ps2a = (h2(i-1,j-1).*csample(1+n-i).*csample(1+n-j));   %quadratic predicted sample
           ps2(n,1) = ps2(n,1) + ps2a';
       end
   end
end
       
psample = ps1(M+1:N/2+M)+ps2(M+1+N/2:N+M);      %predicted sample block

e_block = sample(1+N/2:N)-psample;              %windowed error signal to be concatenated

e = horzcat(e,e_block);                         %concatenate the error block
end

energy_block = e.^2;                            %create energy block for delay minimization
[min_E delay] = min(sum(energy_block(:,:)));    %sum and minimize energy plot and record which delay vector is minimal
tdelay = vertcat(tdelay,delay+10);              %(opional) used to view which delays were sampled for the error signal
error = vertcat(error,e(:,delay));              %construct the error signal using smallest energy vector
Eplot = vertcat(Eplot,energy_block(:,delay));   %construct energy plot using smallest energy vector
countdown = reps-c
end

toc

%generate total Gain of signal given by Energy(s[n])/Energy(e[n])
Energy_e = sum(Eplot);
Gain = sum(Energy_s/Energy_e)              
Gaindb = 10*log10(Gain)