/* gainphase.mac -- Compute the gain and phase offset
* of a sinusoidal function in the time domain or a transfer function in the
* Laplace domain, as a function of frequency and any constants.
*
* copyright 2009 David Simcha
* Released under the terms of the GNU General Public License
*
* Assumptions:
*
* 1. As t -> infinity, the signal approaches a sinusoidal steady state.
* 2. You want the gain and phase as t -> infinity, not at some other arbitrary
* time.
* 3. If you're using gainphase(), you have a full symbolic representation
* of the signal in terms of a frequency variable w and a time variable t.
* If you're working in the Laplace domain, you have a transfer function.
* Examples
(%i1) gainphase(sin(w*t) + cos(w*t) + exp(-t), t, w);
(%o1) [sqrt(2),%pi/4]
(%i2) gainphaselaplace(1 / (s + a), s, w);
(%o2) [1/sqrt(w^2+a^2),-atan2(w/(w^2+a^2),a/(w^2+a^2))]
*/
gainphase(f, t, w) := block(
assume(w >= 0),
gainphaselaplace(
laplace(f, t, s) / laplace(sin(w*t), t, s)
, s, w));
gainphaselaplace(laplaceFunc, s, w) := block([transferFunc, gain, phase, ratprint],
ratprint : false,
assume(w >= 0),
transferFunc : ratsimp(laplaceFunc),
transferFunc : ratsimp(subst(%i*w, s, transferFunc)),
gain : radcan( cabs(transferFunc)),
realTransfer : ratsimp(realpart(transferFunc)),
imagTransfer : imagpart(transferFunc),
phase : if is(realTransfer = 0) then block(
radcan( atan2( imagTransfer, realTransfer)))
else block(
quot : ratsimp( imagTransfer / realTransfer),
atan2( ratnumer(quot), denom(quot))),
[gain, phase]);