%-----------------------------------------------------------------
% Speech Recognition - Numbers 1 to 5
% Energy normalization and time alignment
% References:
% [1] L. Rabiner and B.H. Juang, Fundamentals of Speech
%     Recognition, Prentice-Hall, 1993.
% [2] P.E. Papamichalis, Practical Approaches to Speech
%     Coding, Prentice-Hall, 1987.
% [3] J.D. Markel and A.H. Gray, Linear Prediction of
%     Speech
%-----------------------------------------------------------------
close all
clear all
word = ['one  ';'two  ';'three';'four ';'five '];
swen = 80       % size of window to compute energy
szspeec = 4800  % size of each speech (= of each number)
normenergy = 1; % normalized energy of each speech
swzc = 80       % size window zero crossings
levelzerocross = 0.002  % level zero crossings
swar = 240      % size window autoregressive
load ricardo.mat
load daniel.mat
load mark.mat
for i=1:5       % read training set
   sp(i, 1:10,1:szspeec)=r(i,1:10,1:szspeec);
   sp(i,11:20,1:szspeec)=d(i,1:10,1:szspeec);
   sp(i,21:30,1:szspeec)=m(i,1:10,1:szspeec);
end
for i=1:5       % low-pass signal and
   for j=1:30   % preemphasize speech signal (Ref. [1] pag. 112)
     sp(i,j,:)=sp(i,j,:)-mean(sp(i,j,:));
     sp(i,j,:)=filter([1 -15/16],1,sp(i,j,:));
  end
end
for i=1:5  % plot energy of data set
   figure(1)
   subplot(5,1,i)
   avgg=0;
   for j=1:30
      sp(i,j,:)=sp(i,j,:)./sqrt(sum(sum(sp(i,j,:).^2)))*normenergy;
      for k=1:szspeec/swen
         energy(k)=sum(sum((sp(i,j,1+(k-1)*swen:k*swen).^2)));
      end
      gvct=0;
      for k=1:szspeec/swen
         gvct=gvct+k*energy(k);  % gravity center of gravity
      end
      gvct=gvct/(szspeec/swen)/mean(energy(1:szspeec/swen));
      gvct=round(gvct);
      en=zeros(1,swen);
      tr=40;
      for k=gvct:-1:1  % adust a uniform center of gravity
         en(tr)=energy(k);
         tr=tr-1;
      end
      tr=40;
      for k=gvct:szspeec/swen
         en(tr)=energy(k);
         tr=tr+1;
      end
      figure(1)
      plot(energy,'.')
      axis([1 swen 0 .25]);
      hold on
      figure(2)
      subplot(5,1,i)
      plot(en,'.')
      axis([1 swen 0 .25]);
      hold on
      energy=zeros(1,60);
      vv(i,j)=gvct;  % center of gravity
      s=num2str(j);
      s=strcat(word(i,1:5),'[',s,']');
      disp(s)
   end
   disp(' ')
   ylabel(word(i,1:5));
   figure(1)
   ylabel(word(i,1:5));
end
newsp=[];  % training set with uniform center of gravity
for i=1:5
   for j=1:30
      newsp(i,j,1:szspeec)=zeros(1,szspeec);
      tr=2400;
      for k=vv(i,j)*swen:-1:1
         if(tr>0)
            newsp(i,j,tr)=sp(i,j,k);
            tr=tr-1;
         end         
      end
      tr=2400;
      for k=vv(i,j)*swen:szspeec
         if(tr<=szspeec)
            newsp(i,j,tr)=sp(i,j,k);
            tr=tr+1;
         end         
      end
      s=num2str(j);
      s=strcat('new-',word(i,1:5),'[',s,']');
      disp(s)
      vector=reshape(newsp(i,j,:),1,4800);
   end
   disp(' ')
end
save speech newsp