--- a +++ b/envelop_hilbert.m @@ -0,0 +1,100 @@ +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 \ No newline at end of file