function freq = cutoff(a,b,g)

if ~exist('g','var')
    g = 1/2;
end

o = optimset('display','iter');

freq = lsqnonlin(@objfunc,pi/2,0,pi,o);

    function x = objfunc(p)
        f = frf(a,b,p);
        x = 1000*(abs(f) - g);
    end

end