Fix discrete-time frequency handling in describing_function_response - #1270

For a discrete-time system, describing_function_response found limit
cycles with about the right amplitude but reported the wrong frequency:
it took the oscillation frequency as the imaginary part of the Nyquist
contour, which for a discrete system lies on the unit circle
exp(1j*w*dt), so it reported about sin(w*dt) - the same number whatever
dt is.  The refinement cost function and the returned positions also
evaluated H at z = 1j*w rather than on the unit circle, and refinement
typically failed with "not able to refine result; returning estimate".

For example, H(z) = 0.5/(z^2 - 0.5z) with dt = 1 crosses the negative
real axis at w = 1.3181, but the reported frequency was 0.9703, and
dt = 0.5 gave the same 0.9703 instead of 2.6362.

When H.isdtime(strict=True) - matching the timebase semantics used by
nyquist_response - take the frequencies as np.angle(contour) / dt and
evaluate H at exp(1j*w*dt) in the refinement cost function and for the
returned positions.  The strict check keeps dt=None (unspecified
timebase) on the continuous-time path; dt=True (discrete, unspecified
sample time) reports rad/sample.

Adds regression tests: a known unit-circle crossing for dt = 1 and
dt = 0.5, plus dt=None continuous-equivalence and dt=True coverage.
The discrete tests fail at baseline; continuous-time behavior is
unchanged.

Fixes python-control#1269