% This function filters ventricular peaks by relative amplitude and calculates the amplitude of each peak.

% INPUT PARAMETERS (Input arguments):
% - vec:    Time-series signal containing cardiac events.
% - locs:   Times (in seconds) of the detected peaks.
% - pks:    Amplitude values of the signal at the positions specified by 'locs'.
% - inCua2: Search window around the peak (in number of samples).
% - FreVi:  Sampling frequency of the signal or video (Hz).
% - uhmAm:  Minimum amplitude value a peak can have. It is evaluated with respect to the signal normalized between the minimum and maximum of the detected peaks.
% - tiempo: Time threshold; excludes any peak detected before this value (seconds).

% OUTPUT PARAMETERS (Output arguments):
% - LocAmp: Indices of the peaks that successfully passed the filter.
% - vecAmp: Matrix of dimensions [N x 4] with the structure [AmpI, AmpD, idxIzq, idxDer] for each accepted peak, where:
%   * AmpI:   Peak amplitude with respect to its preceding (left) minimum.
%   * AmpD:   Peak amplitude with respect to its succeeding (right) minimum.
%   * idxIzq: Index of the preceding minimum value for each peak.
%   * idxDer: Index of the succeeding minimum value for each peak.

function [LocAmp,vecAmp]=filterPeaks(vec,locs,pks,inCua2,FreVi,uhmAm,tiempo)
    N=length(vec); % Number of video frames
    LocAmp=[];
    vecAmp=[];
    if tiempo>0 && tiempo<(N/FreVi)
        we=find(locs<=tiempo); %The first second of the signals is eliminated
        inmx=max(we);
    else
        inmx=0;
    end
    t=(1/FreVi)*(1:1:N); %time vector
    hloc=length(locs);
    % Reduced window (80%) for left/right search
    inCua3=round(inCua2*.8);
    VecIzq=[];
    VecDer=[];    
    % Find local minima to the left and right of each valid peak
    for i4=inmx+1:hloc-1
        pos=find(locs(i4)==t);
        if i4==1
            Fin=pos-inCua3;
        else
            Fin=find(locs(i4-1)==t)+round(inCua2*.4);
        end
        AMIzq=min(vec(Fin:pos));
        F1=find(locs(i4+1)==t)-round(inCua2*.4);
        AMDer=min(vec(pos:F1));
    
        posDer=find(vec(pos:F1)==AMDer);
        posIzq=find(vec(Fin:pos)==AMIzq);
        hvi=length(vec(Fin:pos));
        VecIzq=[VecIzq;(pos-(hvi-posIzq(1)))];% Amplitude of the peak to the left
        VecDer=[VecDer;pos+posDer(1)]; % Amplitude of the peak to the Right
    end
    
    % Normalizes using the segment's dynamic range (excludes the first second)
    minV=min(vec(FreVi+1:end));
    maxV=max(vec(FreVi+1:end));
    dmax=maxV-minV; 
    
    % Calculate amplitudes and filter by relative threshold uhmAm
    for i5=inmx+1:hloc-1
        i6=i5-(inmx);
        AmpI=pks(i5)-vec(VecIzq(i6));
        AmpD=pks(i5)-vec(VecDer(i6));
        aref=(AmpI/dmax);
        if aref>=uhmAm
            LocAmp=[LocAmp;i5];
            vecAmp=[vecAmp;AmpI,AmpD,VecIzq(i6),VecDer(i6)];
        end
    end
end