% find 95% containment BW of half-sine pulse
gamma = 0.95;
ff = -20:0.01:20;
P_num = 2*cos(pi*ff);
P_denom = pi*(1-4*ff.^2);

% at ff = +/-0.5, we need to use L'Hospitals rule, P(+/-0.5) = 0.5;
troubleInd = find(abs(abs(ff)-0.5) < sqrt(eps));
P_num(troubleInd) = 0.5;
P_denom(troubleInd) = 1;

% it's safe to divide now
PP = P_num./P_denom;

Sr = cumsum(abs(PP).^2)/sum(abs(PP).^2);

BgInd = find(Sr >= (1+gamma)/2, 1);
Bgamma = 2* ff(BgInd)