function out = VAD_find(x, LWIN, SHIFT, NFFT, FI, NE, THR)

% out = VAD_find(x, LWIN, SHIFT, NFFT, FI, NE)
% 
% This function returns the samples of VAD function of the input signal.
% 
% x         - input signal
% LWIN      - length of win in samples (default: 256)
% SHIFT     - shift in the number of samples (default: 80)
% NFFT      - number of samples used during the FFT calculation 
%             (default: 256)
% FI        - number of frames for noise averaging initizalization
%             (default: 2)
% NE        - noise estimation (default: 0)
% THR       - threshold (default: 1)
% out       - output VAD function

%% Paths and variables
if (nargin < 2 || (isempty(LWIN)))
    LWIN = 256;
end
if (nargin < 3 || (isempty(SHIFT)))
    SHIFT = 80;
end
if (nargin < 4 || (isempty(NFFT)))
    NFFT = 256;
end
if (nargin < 5 || (isempty(FI)))
    FI = 2;
end
if (nargin < 6 || (isempty(NE)))
    NE = 0;
end
if (nargin < 7 || (isempty(THR)))
    THR = 1;
end

Nframes = floor((length(x) - LWIN)/SHIFT) + 1;
SP   = spectrum(x, LWIN, SHIFT, NFFT);
V    = zeros(Nframes, 1);
alfa = 0.98;

%% Apply VAD
T = zeros(Nframes, 1);

for frame = 1:Nframes       
    if (frame <= FI)
       NE = (1 - 1/frame)*NE + 1/frame*SP(:, frame);
       V(frame) = 0;
    else
       T(frame) = 10*log10(1/(NFFT/2)*sum(SP(:, frame).^2./NE.^2));
       if (T(frame) > THR)
           V(frame) = 1; 
       else
           V(frame) = 0;
           NE = alfa*NE + (1 - alfa)*SP(:, frame);
       end
    end
end

%% Calculate the output
y    = x;
out  = x;
xmax = max(x);

N = min(floor((length(x) - LWIN)/SHIFT) + 1, length(V));

for i = 1:N
	if (V(i) == 0)
        y((i - 1)*SHIFT + [1:LWIN])   = 0;
        out((i - 1)*SHIFT + [1:LWIN]) = 0;
    else
        y((i - 1)*SHIFT + [1:LWIN])   = 1.1*xmax;
        out((i - 1)*SHIFT + [1:LWIN]) = 1;
	end
end

%% Plot signal and VAD if necessary
if (nargout == 0)
    plot([x y]);
    grid on;
    H = axis;
    axis([H(1) length(x) H(3) H(4)]);
end

function SP = spectrum(x, LWIN, SHIFT, NFFT)

% This function calculates the spectrum of each frame in input signal.
% 
% x         - input signal
% LWIN      - length of win in samples
% SHIFT     - shift in the number of samples
% NFFT      - number of samples used during the FFT calculation
% SP        - the output spectrum

L = floor((length(x) - LWIN)/SHIFT) + 1;
SP = zeros(NFFT, L);

for i = 1:L
    r = (i - 1)*SHIFT + [1:LWIN];
    s_w = x(r);
    SP(:, i) = abs(fft(s_w, NFFT));
end