function J = integrate_freqresp(sys, f)
    arguments
        sys (1,:) lti
        f (1,2) double {mustBePositive} = [1e0,1e6];
    end
    
    J = sqrt(integral(@(f) squeeze(sum(abs(freqresp(sys, f, "Hz")).^2, 2))', f(1), f(2)));
end

