%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extrapolation Algorithm - Minimum Weighted Norm Extrapolation   %

% with Frequency Domain Blocking,                                  %

% RMS factor, and cross-fade overlapping                           %

%                                                                  %

% Created by J. Alex Souppa (4/4/99)                               %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

% Algorithm variables

L = 2048;                       % Length of known portion of the signal

B = 256;                        % Size of each block

Overlap = 1/16;                 % Block overlap percentage; 1/16=6.25%

% Computed constants

N = 2*L;                        % Size of input: (known + extrapolated)

Bdiv2 = B/2;                    % Half of block

Onum = Bdiv2*Overlap;          % Number of samples of overlap

K = ceil(L/(Bdiv2-Onum));      % Number of blocks

% Cross-fade Windows

Crossfade=triang(Onum*2)';

CFW(1:Bdiv2)=ones(size(Bdiv2));

CFW(Bdiv2-Onum+1:Bdiv2)=Crossfade(Onum+1:2*Onum);

CFW1=CFW;

CFW(1:Onum)=Crossfade(1:Onum);

CFW(Bdiv2+1:B)=CFW(1:Bdiv2);

CFW1(Bdiv2+1:B)=fliplr(CFW1(1:Bdiv2));

%%% Sine wave test signal

Fs = 8000;                      % Sampling frequency for sine wave sample

Start_Hole = L+1;               % Starting sample of dropout

t=0:1/Fs:((N-1)/Fs);           % Time function

f=sin(2*pi*350*t)+0.5*sin(2*pi*1000*t)+0.7*sin(2*pi*2000*t)+

0.2*sin(2*pi*3000*t);

% Fill the buffer with known data

g=f(Start_Hole-L:Start_Hole-1);

h(1:L)=g;

h(N)=0;                            % Zero-pad the input

% Setup timer

start_time = get_time;

% Frequency Domain Blocking

Fhz=fft(h,N);

dc=Fhz(1);

Nyq=Fhz(L+1);

for x=1:K

   % Block divide input spectrum

   FhBlockz(x,1:Bdiv2)=Fhz((Bdiv2-Onum)*(x-1)+1:(Bdiv2*x)-Onum*(x-1));

   % Zero-pad the last block

   if (x==K)

      temp=(Bdiv2*x)-Onum*(x-1)-L;

      if (temp~=0)

         FhBlockz(x,Bdiv2-temp+1:Bdiv2)=zeros(size(temp));

      end

   end

   % Fill in second half of block

   FhBlockz(x,Bdiv2+2:B)=conj(fliplr(FhBlockz(x,2:Bdiv2)));

   % Define time-domain blocks

   hBlock(x,1:B)=real(ifft(FhBlockz(x,:),B));

   gBlock(x,1:Bdiv2)=hBlock(x,1:Bdiv2);

end

% ******************* Compute extrapolation **********************

% Compute circular autocorrelation

for x=1:K

   o(x,1:B)=cir_auto(hBlock(x,:));

end

% Form matrix G

AA=[];

for ii=1:Bdiv2

   AA=[ii AA];

   for x=1:K

      G(ii,1:Bdiv2,x)=o(x,[AA B:-1:B-Bdiv2+ii+1]);

   end

end

% Compute the vector of constants

for x=1:K

   tempG(1:Bdiv2,1:Bdiv2)=G(:,:,x);

   tempg=gBlock(x,:);

   if (sum(tempg)==0)

      a(x,1:Bdiv2)=ones(size(Bdiv2));

   else

      a(x,1:Bdiv2)=((tempG-o(x,1)*10^(-8)*eye(Bdiv2))\tempg')';

   end

end

% Form the estimate

for x=1:K

   fBlock(x,1:Bdiv2)=gBlock(x,:);

end

for ii=Bdiv2+1:B

   for x=1:K

      fBlock(x,ii)=sum(a(x,:).*o(x,[ii:-1:ii-Bdiv2+1]));

   end

end

% Adjust the signal with RMS factor

for x=1:K

   % RMS of known signal

   RMS_known=sqrt(sum(fBlock(x,1:Bdiv2).^2))/Bdiv2;

   % RMS of extrapolation

   RMS_extrap=sqrt(sum(fBlock(x,Bdiv2+1:B).^2))/Bdiv2;

   % Form output after AutoGain

   fBlockA(x,1:Bdiv2)=fBlock(x,1:Bdiv2);

   fBlockA(x,Bdiv2+1:B)=(RMS_known/RMS_extrap)*fBlock(x,Bdiv2+1:B);

end

% Inverse Frequency Domain Blocking

FfEst=zeros(size(1:N));

for x=1:K

   % Get spectrum of time-domain block

   FfBlockA(x,1:B)=fft(fBlockA(x,:),B);

   % Zero-pad last block

   if (x==K)

      temp=(Bdiv2*x)-Onum*(x-1)-L;

      if (temp~=0)

         FfBlockA(x,Bdiv2-temp+1:Bdiv2)=zeros(size(temp));

      end

   end

   % Window the spectrum block

   if (x==1)

      FfBlockA(x,1:B)=FfBlockA(x,1:B).*CFW1;

   else

      FfBlockA(x,1:B)=FfBlockA(x,1:B).*CFW;

   end

   % Sum overlapping blocks

   if (x~=1)

      FfBlockA(x,1:Onum)=FfBlockA(x,1:Onum)+FfEst((Bdiv2*(x-1))- Onum*(x-1)+1:(Bdiv2*(x-1))-Onum*(x-2));

      FfBlockA(x,B-Onum+2:B)=conj(fliplr(FfBlockA(x,2:Onum)));

   end

   % Replace block in original position of main spectrum

   if (x==K)

      FfEst((Bdiv2-Onum)*(x-1)+1:L)=FfBlockA(x,1:Bdiv2-temp);

   else

      FfEst((Bdiv2-Onum)*(x-1)+1:(Bdiv2*x)-Onum*(x1))= FfBlockA(x,1:Bdiv2);

   end

end

% Fill in second half of main spectrum

FfEst(L+2:N)=conj(fliplr(FfEst(2:L)));

FfEst(1)=dc;

FfEst(L+1)=Nyq;

% Get time-domain equivalent = main output of system

fEst=real(ifft(FfEst,N));

% Replace first half with the known data

fEst(1:L)=g;

% Compute computational time of extrapolation

time=get_time-start_time;

% Find error signal

E=f(Start_Hole:Start_Hole+L-1)-fEst(L+1:N);

E2=abs(E).^2;

MSE=sum(E2)/L;

% Display the results

disp('*********************Results**************************')

if (N==B)

   disp(sprintf('FFT size=%d No Blocking',N))

else

   disp(sprintf('FFT size=%d Block Size=%d (with %.2f%% overlap)',N,B,Overlap*100))

end

disp(sprintf('Computational Time = %.2f seconds',time))

disp(sprintf('Mean Square Error = %f',MSE))