%
%                       PROGRAM SPECT
%                (c) 2000. Rodrigo Quian Quiroga
%
% Given a 1 column input ascii file, calculates its Windowed Fourier 
% Transform and then calculates the power spectra and  the Shannon 
% and Kullback-Leibler (relative) entropies. It also calculates the 
% relative band ratio (RIR).
%
% Data file name, sampling rate, maximum frequency to be considered 
% for the calculus of the entropy, window size, window overlapp and 
% reference window should be given.
% Axis ranges should be changed according with the data.
% The Kullback-Leibler (variable kl) and Shannon (variable h) entropies 
% can be saved with a save command (see matlab help).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
x = load('130-2_c4.txt');   %input file
lx = length(x);
sr = 102.4;              %sampling rate in Hz
ws = 512;              %data window size
ov = 0.5;               %windows overlapp
ref = 1;                %reference window for the KL entropy (in sec)
smooth = 3;             %half size of the spectral smoothing window
fmax =50;               %maximum frequency considered for the entropies

% REFERENCE SPECTRUM
clear xref
  for i=1:ws;
   xref(i) = x(i+floor(ref*sr));
  end
   Y = fft(xref);
   Y(1)=[];
  for i=1:fmax*ws/sr + smooth;
   powerref(i) = abs(Y(i)).^2;
  end
 %frequency window
   totalpower = 0;
    for i=1:fmax*ws/sr;
      powerrefw(i) = 0; 
       for j=-smooth:smooth;         %Barlett-Priestley freq. window  
        if j>-i;
         powerrefw(i) = powerrefw(i) + powerref(i+j) * (1-(j/smooth)^2);
        end
       end
      totalpower = totalpower + powerrefw(i);
    end
    for i=1:fmax*ws/sr;
      powerrefw(i) = powerrefw(i)/totalpower;
    end


%
% LOOP OVER WINDOWS
%

for k=1:floor(lx/(ws*ov))-1;
 clear xt
  for i=1:ws;
   xt(i) = x(i+(k-1)*ws*ov);
  end

% FOURIER AND POWER
Y = fft(xt);
Y(1)=[];
for i=1:fmax*ws/sr + smooth;
 power(k,i) = abs(Y(i)).^2;
end

%BARLETT-PRIESTLEY FREQUENCY WINDOW
totalpower = 0;
 for i=1:fmax*ws/sr;
   powerw(k,i) = 0; 
    for j=-smooth:smooth;
      if j>-i;
       powerw(k,i) = powerw(k,i) + power(k,i+j) * (1-(j/smooth)^2);
      end
    end
   totalpower = totalpower + powerw(k,i);
 end

% ENTROPY 
h(k) = 0; kl(k) = 0;
delta(k) =0;theta(k)=0;alfa(k)=0;beta(k)=0;gamma(k)=0;
for i=1:fmax*ws/sr;
 powerw(k,i) = powerw(k,i)/totalpower;
  if powerw(k,i) > 0; 
   h(k) = h(k) - powerw(k,i) * log2(powerw(k,i));
    if powerrefw(i) > 0;
     kl(k) = kl(k) - powerrefw(i)*log2(powerrefw(i)/powerw(k,i));
    end
  end
  if i< 3.5 * ws/sr; delta(k) = delta(k) + powerw(k,i); 
   elseif i< 7.5 * ws/sr; theta(k) = theta(k) + powerw(k,i); 
    elseif i< 12.5 * ws/sr; alfa(k) = alfa(k) + powerw(k,i); 
     elseif i< 30 * ws/sr; beta(k) = beta(k) + powerw(k,i); 
    elseif i< 60 * ws/sr; gamma(k) = gamma(k) + powerw(k,i); 
  end


end
end

data = (0:ws*ov/sr:(k-1)*ws*ov/sr);
freq = ((1:fmax*ws/sr)/(fmax*ws/sr)) * fmax;
%load datasignal
figure
subplot(5,1,1) 
plot(x) 
axis([0 length(x) -2000 2000]);grid; 
ylabel('x_n')
subplot(5,1,3)
plot(data,h)
axis([0 lx/sr 0 10]);grid; 
ylabel('H')
subplot(5,1,4)
plot(data,kl)
axis([0 lx/sr -10 0]);grid; 
ylabel('K-L')
subplot(5,1,2)
contour(data,freq,powerw',16)
ylabel('f (Hz)')
axis([0 lx/sr 0 fmax]);grid;
subplot(5,1,5)
plot(data,delta,data,theta,data,alfa,data,beta,data,gamma)
axis([0 lx/sr 0 1]);grid; 
ylabel('RIR')