/* get important frequency values for continuous system - see freqs.spl */
freqsvals(a, b, npts)
{
        local i, integ, w, xi, f, z;

        if (argc == 2)
        {
                // assume cascade
                (ep, tz, k) = cas2zp(a);

                npts = b;
        }
        else
        {
                // direct form
                ep = roots(b);
                tz = roots(a);
        }

        if (isempty(ep))
        {
                ep = -1000;
        }

        // does not handle zeros greater than 1e8
        ez = {ep[find(imag(ep) >= 0)], tz[find(abs(tz) < 1e8 && imag(tz) >= 0)]};

        // round start and end to nearest decade
        integ    = abs(ez) < 1e-10; // cater for systems with pure integrators
        highfreq = round(log10(max(3 * abs(real(ez) + integ) + 1.5 * imag(ez))) + 0.5);
        lowfreq  = round(log10(0.1 * min(abs(real(ez + integ)) + 2 * imag(ez))) - 0.5);

        // base frequency range
        diffzp = length(ep) - length(tz);
        w      = logspace(lowfreq, highfreq, npts + diffzp + 10 * (sum(abs(imag(tz)) < abs(real(tz))) > 0));
        ez     = ez[find(imag(ez) > abs(real(ez)))];

        // poles and zeros for oscillations
        if (not(isempty(ez)))
        {
                f     = w;
                npts2 = 2 + 8 / ceil(abs((diffzp + eps) / 10));

                (dum, ind) = sort(-abs(real(ez)), 1);

                z = {};

                loop (i = ind)
                {
                        r1 = max({0.8*imag(ez[i]) - 3 * abs(real(ez[i])), 10 ^ lowfreq});
                        r2 = 1.2 * imag(ez[i]) + 4 * abs(real(ez[i]));
                        z  = z[find(z > r2 || z < r1)];
                        f  = f[find(f > r2 || f < r1)];
                        z  = {z, logspace(log10(r1), log10(r2), sum(w <= r2 && w >= r1) + npts2)};
                }

                w = sort({f, z}, 1);
        }

        /* interpolate on a uniform log grid */
        xi = linspace(1, length(w), npts - 1);
        w  = xylookup(1..length(w), log10(w), xi);

        /* add DC */
        w  = {0, 10^yvals(w)};

        return(w);
}