function res = identify_signal_features(x, fs, thr_level_mV)
% %% Function identifies features of AE signal for futher analysis. 
% This script uses a ?oating threshold in order to remove part of the signal
% due to reflection at specimen edges. 
%
% [AMP,DUR,RIS,CNT,CTM,KUR,ENG,FR1,FR2,CDD,FAM,FCG]  = identify_ae_features(x, fs, level)
%
%
%          x   - input signal vector
%          fs  - sampling frequency (default = 10 MHz)
% thr_level_mV - percentage threshold for determination several features 
%     (default = 0.2 <= 20%)
%
% td_amp_max    - maximalni amplituda
% td_dur        - DURATION
% td_ris        - RISE time
% td_rise_angle - nastupni uhel udalosti
% td_amp_first  - velikost amplituda prvniho prekmitu
% td_cnt        - COUNT
% td_ctm        - COUNT to PEAK
% td_zcr        - pocet pruchodu nulou
% td_cdd        - pomer hodnot pocet pruchodu nulou/ trvani AE udalosti
% td_eng        - energie signalu nad urovni
% td_kur        - statisticky moment 4 radu
% td_skw        - statisticky moment 3 radu
% 
% fd_amp_max    - maximalni amplituda ve frekvencni oblasti
% fd_freq_max   - frekvence odpovidajici maximu ve frekvencni oblasti
% fd_centrumgravity - frekvence teziste FFT signalu
% fd_power1     - energie pro pasmo  10-110kHz ve frekvencni oblasti
% fd_power2
% fd_power3
% fd_power4
% fd_power5
% fd_power6 
% 
% tf_freq_first - frekvence prvniho peaku ve CWT
% tf_freq_max   - frekvence maximalniho peaku v CWT
% 
%  ------------------------------------------------------------------------------
%  Petr Sedlak, 2010, Department of physics, Brno University of Technology

% check input
if (nargin < 1)|(length(x) == 1), error('The input have to be a vector.');end
if (nargin < 2), fs = 1e7; end
if (nargin < 3), thr_level_mV = 10; end
if (thr_level_mV >= 1000)|(thr_level_mV < 10)
    thr_level_mV = 10; 
    disp('rozhodovaci uroven by mela byt v intervalu 10 - 1000 mV');
end

td_amp_max=0;td_dur=0;td_ris=0;td_ris_angle=0;td_amp_first=0;td_cnt=0;td_ctm=0;td_zcr=0;td_cdd=0;td_eng=0;
td_kur=0;td_skw=0;fd_amp_max=0;fd_freq_max=0;fd_centrumgravity=0;fd_power1=0;fd_power2=0;fd_power3=0;
fd_power4=0;fd_power5=0;fd_power6=0;tf_freq_first=0;tf_freq_max=0;
flg_noise = 0; % priznak sumove udalosti

% urci zakladni parametry 
N = length(x);      % delka signalu
t = [0:(N-1)]./fs;  % vektor casu
W  = tukeywin(N,0.1); % vytvor okno

DELAY = 400; % difference in sampling points which is considered as the end of event

%% -- PARAMETRY V CASOVE OBLASTI ------------------------------------------

% maximalni amplituda
td_amp_max = max(abs(x)); % % vysvetleni zkratky time domain% tf_freq_first - frekvence prvniho peaku ve CWT
% tf_freq_max   - frekvence maximalniho peaku v CWT - amplitude - maximal


% urci trvani udalosti
idxes_event = find(abs(x) >= thr_level_mV);
AA = find(diff(idxes_event)> DELAY,1);
if isempty(AA), AA = length(idxes_event); end
idxes_event = idxes_event(1:AA);
if ~isempty(idxes_event)
   td_dur = (idxes_event(end)-idxes_event(1))/fs;
else
    flg_noise = 1;
end

if flg_noise  == 0
    % urci rise time
    [a,idx_MAX] = max(x);
    td_ris = (idx_MAX-idxes_event(1))/fs;
    
    % urci rise angle
    td_ris_angle = atan(td_amp_max/(td_ris*1e6))*180/pi; % pocitej cas v mikrosekundach
    
    % urci amplitudy prvni spicky signalu
    y = abs(x((idxes_event(1)-1):idxes_event(end)));
    idx = find(diff(sign(diff(y)))==-2,1) +  idxes_event(1) - 1;
    td_amp_first = abs(x(idx));
    
    % urci pocet period prochazejici urovni (COUNT)
    y = x -  thr_level_mV;
    idxes_CNT = find(diff(sign(y((idxes_event(1)-1):idxes_event(end)))) == 2)+idxes_event(1)-1;
    td_cnt = length(idxes_CNT);  % counts
    
    % urci counts to peak
    idxes_CTM = idxes_CNT(idxes_CNT <= idx_MAX);
    td_ctm = length(idxes_CTM);  % counts to peak
    
    % urci pocet zero crossing ZCR
    idxes_ZCR = find((abs(diff(sign(x((idxes_event(1)-1):idxes_event(end))))) == 2))+idxes_event(1)-1;
    td_zcr = length(idxes_ZCR)/2;
    
    % urci counts/duration
    td_cdd = td_zcr/td_dur;
    
    % urci energii
    env = abs(calculate_envelope(x(:).*W(:)));% konstrukce obalky
    idx_ENG = find(env >= thr_level_mV);
    env = env - thr_level_mV;
    td_eng = sum(env(idx_ENG));
end

% urci kurtosis
td_kur = sum(((x-mean(x))/std(x)).^4)/N - 3;

% urci skewness
td_skw = sum(((x-mean(x))/std(x)).^3)/N;

%% -- PARAMETRY VE FREKVENCNI OBLASTI -------------------------------------
% vypocet ve frekvencni oblasti ma zakladni predpoklad, ze vzorkovaci
% frekvence byla 10MHz, a proto muzeme vypocet zprumerovat 4x
% respektive snizi vzorkovaci frekvenci 4x
NN = 4;
NFFT = 2^(nextpow2(N/NN));
L = round(length(x)/NN);
Y = zeros(NFFT/2,1);
for k = 1:NN
    y = x(k:NN:end);
    y_fft = fft(y,NFFT)./N;
    Y = Y +  2*(abs(y_fft(1:NFFT/2)))';
    f = fs/NN /2*linspace(0,1,NFFT/2);
end
Y = Y./NN;

% urci maximalni amplitudu ve frekvencni oblasti
[a,i] = max(Y);
fd_amp_max = a;

% urci frekvenci odpovidajicic maximu
fd_freq_max = f(i);

% urci spectrum's center of gravity
fd_centrumgravity = sum(Y(:).*f(:))./sum(Y);

% urci energii pro pasmo  10-110kHz
ii = find( (f >= 10e3)&(f < 110e3));
fd_power1 = sum(Y(ii))/sum(Y);

% urci energii pro pasmo 110-210kHz
ii = find( (f >= 110e3)&(f < 210e3));
fd_power2 = sum(Y(ii))/sum(Y);

% urci energii pro pasmo 210-310kHz
ii = find( (f >= 210e3)&(f < 310e3));
fd_power3 = sum(Y(ii))/sum(Y);

% urci energii pro pasmo 310-410kHz
ii = find( (f >= 310e3)&(f < 410e3));
fd_power4 = sum(Y(ii))/sum(Y);

% urci energii pro pasmo 410-510kHz
ii = find( (f >= 410e3)&(f < 510e3));
fd_power5 = sum(Y(ii))/sum(Y);

% urci energii pro pasmo 510-610kHz
ii = find( (f >= 510e3)&(f < 610e3));
fd_power6 = sum(Y(ii))/sum(Y);

%% -- PARAMETRY VE CASOVE-FREKVENCNI OBLASTI ------------------------------
if flg_noise == 0
    
    f = linspace(1e3,8e5,200); % urci frekvenci
    % vypocet vlnkovy transformace
    WT = ContinuousWaveletTransform(x,fs,f,'gabr');
    
    % najdi localni maxima
    Locations = MatrixLocalMax(abs(WT));
    [iif,iit] = find(Locations == 1);
    % for k = 1:length(iif)
    %     K(k) = abs(WT(iif(k),iit(k)));
    % end
    K = diag(abs(WT(iif,iit)));
    [K, ii] = sort(K,'descend');
    
    ii = ii(K > 0.5);
    iif = iif(ii);
    iit = iit(ii);
    
    % urci frekvenci prvniho peaku
    [a,i] = min(iit);
    tf_freq_first = f(iif(i));
    
    % urci frekvenci maximalniho peaku
    tf_freq_max = f(iif(1));
end

res = [td_amp_max, td_dur, td_ris, td_ris_angle, td_amp_first,td_cnt, td_ctm, td_zcr, td_cdd, td_eng, td_kur, td_skw, fd_amp_max, fd_freq_max, fd_centrumgravity, fd_power1, fd_power2,fd_power3, fd_power4, fd_power5, fd_power6,tf_freq_first, tf_freq_max];
return

%% calculate envelope 
function env = calculate_envelope(xr)
% function calculate envelope of the signal on the basis of HIlbert
% transform

if ~isreal(xr)
  xr = real(xr); hilbert
end
% Work along the first nonsingleton dimension
[xr,nshifts] = shiftdim(xr);
n = size(xr,1);

u = fft(xr,n,1); % n-point FFT over columns.
h  = zeros(n,~isempty(u)); % nx1 for nonempty. 0x0 for empty.
if n > 0 && 2*fix(n/2) == n
  % even and nonempty
  h([1 n/2+1]) = 1;
  h(2:n/2) = 2;
elseif n>0
  % odd and nonempty
  h(1) = 1;
  h(2:(n+1)/2) = 2;
end
u = ifft(u.*h(:,ones(1,size(u,2))));

% Convert back to the original shape.
u = shiftdim(u,-nshifts);

env = abs(u);
return

%% find maxima
function Locations = MatrixLocalMax(A,varargin)
%MATRIXLOCALMAX Find local maxima in a matrix by returning a
% new matrix that equal one at the maxima and zero elsewhere.
%
% Locations = MatrixLocalMax(A) finds all isolated local maxima
% of the matrix A.
%
% Locations = MatrixLocalMax(A,Reach) finds local maxima of the
% matrxi A that are separated by at least Reach points. Reach
% should be an integer larger than or equal to one. Reach is
% included to filter maxima in matrices with noise.
%
% Note: If two points within Reach distance of one another are
% identical, MatrixLocalMax only returns one value
%
% Example:
% A = peaks;
% Locations = MatrixLocalMax(A);
%
% See also MAX


%
% Author: David Eyre
% University of Utah
% Copyright 2000
%
if nargin==1
    Reach = 1;
elseif nargin==2
    Reach = max(round(varargin{1}),1);
else
    error('Incorrect number of input arguments')
end


[m,n] = size(A);
P = (1:m)+Reach;
Q = (1:n)+Reach;

%
% Pad A so boundary conditions become transparent
%
APad = (min(min(A))-1)*ones(m+(2*Reach),n+(2*Reach));
APad(P,Q) = A;
   
%
% Line up arrays to determine local maxima in boxes centered
% around every point with radius = Reach. This uses a lot of
% memory but is more efficient that non-vector approaches
%
Count = 0;
Sheets = (2*Reach+1)^2;
Work = zeros(m,n,Sheets);
for j=-Reach:Reach
    for k=-Reach:Reach
        Count = Count+1;
        Work(:,:,Count) = APad(P+j,Q+k);
    end
end
[WorkMax,IMax] = max(Work,[],3);
   
%
% Maxima correspond to points mapped to themselves
%
Invariant = ((2*Reach+1)^2+1)/2;
Locations = 1-sign(abs(IMax-Invariant));
return
%% calculate continuous wavelet transform
function [res] = ContinuousWaveletTransform(x, fsample, frange, wname)
% function calculates continuous wavelet transform in frequency domain
% The possible mother wavelets are :
%   - Morlet 'morl'
%   - Gabor 'gabr' (default)
%   - Mexican Hat 'mexh'
%   - Battle-Lemarie Wavelet 'bale'
%
% INPUT VARIABLES
%       x - time series (signal)
% fsample - sampling frequency of signal
%  frange - frequencies of user's interest
%   wname - name of mother wavelet ('morl','gabr','mexh')
%
% OUTPUT VARIABLES
%     res - 2D time-frequency distribution of CWT coefficients
%
% -----------------------------------------------------------
% Petr Sedlak (c), The University of Tokyo 2009

% check input variables
if nargin < 3, error('Not enough input parameters. :-('); elseif nargin < 4, wname='gabr'; end

% compute the angular frequency omega_k 
n  = length(x); n_2 = round(n/2); 
omega_k = 2*pi.*[((1:n_2)-1)./n -(n-((1+n_2):n)+1)./n];
omega_k(1)=1e-9;

% fast fourier transform of signal
x_fft = fft(x);

% calculate scale
s = fsample./frange;
ns = length(s); 

switch wname % choice of mother wavelet
    case 'morl'
        wavelet_fun = @(w) (pi^-0.25).*exp(-0.5.*((w - 2*pi).^2));
    case 'mexh'
        wavelet_fun = @(w) (-((w).^2).*exp(-(w).^2));
    case 'gabr'
        wavelet_fun = @(w) (sqrt(2*pi)*(pi^-0.25)*sqrt(sqrt(2/log(2))/2)).*exp(-0.5*(0.5/log(2)).*((w-2*pi).^2));
    case 'bale'
        wavelet_fun = @(w) ((1/sqrt(2*pi)).*exp(-sqrt(-1).*w./2).*[(sin(w./2)./(w./2)).^2]);
end

% the algorithm
for is = 1:ns
    z = s(is).*omega_k;
    xg = sqrt(2*pi*s(is)).*wavelet_fun(z)+0*i;
    xcmor(is,:) = ifft(xg(:).*x_fft(:));    
end

% norm the matrix
res =  xcmor./max(max(abs(xcmor)));