function [evm_total,evm_lvl] = calc_evm(test_signal,reference_signal,options)

arguments(Input)
    test_signal;
    reference_signal;
    options.skip_front = 0;
    options.skip_end = 0;
    options.returnErrorLocation = 0;
end

options.skip_end = abs(options.skip_end);
options.skip_front = abs(options.skip_front);

assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");

if isa(reference_signal,'Signal')
    reference_signal = reference_signal.signal;
end

if isa(test_signal,'Signal')
    test_signal = test_signal.signal;
end

% TRIM
[test_signal,reference_signal]=trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);

% CALC EVM
[evm_total,evm_lvl] = calc_evm_(test_signal,reference_signal);


    function [evm_total, evm_lvl] = calc_evm_(test_signal, reference_signal)
        % Validate input
        assert(length(test_signal) == length(reference_signal), "Sequence length does not match");

        % Calculate error vector
        error_vector = test_signal - reference_signal;

        % EVM (RMS) as percentage, per MathWorks definition
        evm_total = sqrt(sum(abs(error_vector).^2) / sum(abs(reference_signal).^2)) * 100;

        % Per-level EVM
        k = unique(reference_signal);
        evm_lvl = NaN(1, length(k));
        for lvl = 1:length(k)
            idx = reference_signal == k(lvl);
            if any(idx)
                lvl_errors = error_vector(idx);
                lvl_refs = reference_signal(idx);
                evm_lvl(lvl) = sqrt(sum(abs(lvl_errors).^2) / sum(abs(lvl_refs).^2)) * 100;
            end
        end
    end

    function [data_,reference_]=trimseq(data,reference,skipstart,skip_end)

        data_  =  data(skipstart+1:end-skip_end,:);

        delta_bits = length(reference) - length(data);

        skip_end = delta_bits + skip_end;

        reference_ = reference(skipstart+1:end-skip_end,:);

    end

end