function alarm = envelop_hilbert(y,Smooth_window,threshold_style,DURATION,gr)

% y = Raw input signal to be analyzed

% Smooth_window :this is the window length used for smoothing your signal in terms of number of samples

% threshold_style : set it 1 to have an adaptive threshold OR set it 0 to manually select the threshold from a plot

% DURATION : Number of the samples that the signal should stay in terms of number of samples

% gr = make it 1 if you want a plot and 0 when you dont want a plot

% Tuning parameters for the best results;

% 1. DURATION is correlated to your sampling frequency, you can use a multiple of your sampling frequency e.g. round(0.050*SamplingFrequency)

% 2. Smooth_window is correlated to your sampling frequency, you can use a multiple of your sampling frequency e.g. round(0.0500*SamplingFrequency), this is the window length used for smoothing your signal

% alarm : vector resembeling the active parts of the signal

%Method

% Calculates the analytical signal with the help of hilbert transfrom,

% takes the envelope and smoothes the signal. Finally , with the help of an

% adaptive threshold detects the activity of the signal where at least a

% minimum number of samples with the length of 

% (DURATION) Samples should stay above the threshold). The threshold is a

% computation of signal noise and activity level which is updated online.

%% Input handling

%--------------------------- input handling ---------------------------%

if nargin < 5

    gr = 1;

    if nargin < 4

        DURATION = 20;                                                     % default

        if nargin < 3

           threshold_style = 1;                                            % default 1 , means it is done automatic

           if nargin < 2

             Smooth_window = 20;                                           % default for smoothing length

             if  nargin < 1

               v = repmat([.1*ones(200,1);ones(100,1)],[10 1]);            % generate true variance profile

               y = sqrt(v).*randn(size(v));

             end

           end

        end

    end

end

%% calculate the analytical signal and get the envelope

test=y(:);

analytic = hilbert(test);

env = abs(analytic);

%% take the moving average of analytical signal

env = conv(env,ones(1,Smooth_window)/Smooth_window);                       % smooth

env = env(:) - mean(env);                                                  % get rid of offset

env = env/max(env);                                                        % normalize

%% threshold the signal

if threshold_style == 0

   hg=figure;plot(env);title('Select a threshold on the graph')

   [~,THR_SIG] =ginput(1);

   close(hg);

end

% ------------------------- Threshold Style ---------------------- %

if threshold_style

   THR_SIG = 4*mean(env);

end

nois = mean(env)*(1/3);                                 % noise level

threshold = mean(env);                                  % signal level

% ------------------- Initialize Buffers -------------------------%

thres_buf  = zeros(1,length(env)-DURATION);

nois_buf = zeros(1,length(env)-DURATION);

THR_buf = zeros(1,length(env));

h=1;

alarm =zeros(1,length(env));

for i = 1:DURATION/2:length(env)-DURATION

  %------ update threshold 10% of the maximum peaks found ------------%

  if env(i:i+DURATION) > THR_SIG 

      alarm(i) = max(env);

      threshold = 0.1*mean(env(i:i+DURATION)); 

      h=h+1;

  else

      % ----------- update noise --------------%

      if mean(env(i:i+DURATION)) < THR_SIG

          nois = mean(env(i:i+DURATION)); 

      else

          if ~isempty(nois_buf)

              nois = mean(nois_buf);

          end

      end

  end 

  thres_buf(i) = threshold;

  nois_buf(i) = nois;

  % ------------------------- update threshold -------------------------- %

  if h > 1

     THR_SIG = nois + 0.50*(abs(threshold - nois));            

  end

  THR_buf(i) = THR_SIG;

end