Fix discrete-time frequency handling in describing_function_response - #1272

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 #1269

written by claude opus 5.5, reviewed by me