Algorithmically get Amplitude and Phase of Sine wave?

Question

I'm trying to figure out a way to algorithmically get the amplitude and phase of a function that has sinusoidal terms in the Maxima computer algebra system. This only applies to steady state (as t -> infinity and the transients decay). For example, a trivial case would be:

f(t) = 1 / w * sin(w * t + theta) + exp(-a * t) + 8

In this case, the gain would be 1 / w, the phase offset would be theta, and we would ignore the transient term exp(-a * t) because we only care about the steady state gain and phase delay, and exp(-a * t) -> 0 as t -> infinity. We would also ignore the "+ 8" term because it's just a DC offset. The way I've been taught to do it in my engineering classes requires a lot of heuristics and tedious rearranging of equations to get them in a form similar to the above, where the answer is obvious just from looking at it.

Does anyone know of a general, algorithmic method of finding the gain and phase delay assuming they exist, given that I have the full power of a computer algebra system (and the standard functions that one would expect a CAS to have) to throw at it? Although I will likely be implementing it in Maxima, I would certainly appreciate generic answers explained just in terms of math.

Edit: I thought it was abundantly clear from my example that I want the answer symbolically, in terms of w. w is really supposed to be omega, and represents the frequency of the input. What I'm really asking is whether there are any standard math operations that will produce the gain and phase terms without a bunch of heuristic, manual equation rearranging.

Answer 1

Do you mean symbolically or numerically?

Numerically you'd want to perform a Fourier transform:

• sample the function for at at least twice the expected maximum frequency (even higher if you want more precise measurement of phase) and for as long as at least your maximum expected wavelength

• perform a Fourier transfomr (searching for FFT should turn up lots of examples - my searches suggest that maxima may even have a built in fft function)

• this'll give you a "frequency domain" version of your function. You'll have a series of complex values where the magnitude is the amplitude for that frequency, and the angle is the phase of that frequency component. In your case it sounds like you want to look for the frequeny with the peak amplitude

Answer 2

This may or may not be ideal, but I'm assuming you can't perform any math on the function generating the waveform (or the function is not available) -

I would sample the wave at intervals that are significantly less than the period. This is probably difficult - you'll need to know for certain that the wave's peirod is greater than the sampling interval, but you won't want it to be too small or it will take forever to finish a few cycles of the wave. Perhaps you could start sampling with some minimal value, and ramp it up until the value changes at some acceptable rate.

If you do this, then you can shift the wave to center on the X-axis (by subtracting or adding the average value) and determine 1) Min/max values, which provides an amplitude and 2) X-intercepts, which provides the period.

Answer 3

Either a Laplace or a Fourier transform may help; you can apply either symbolically as well as numerically. But I'm not sure you will be able to create a general algorithmic method - usually there are different case depending where the poles are.

Answer 4

After thinking about this for a while, I've figured this out myself. I'll post it here because it's pretty interesting. The assumptions are:

1. The signal has existed for infinite time, so everything is in steady state.
2. Your expression only has one frequency in it. This is easy to check just by inspection, even if it is in some horribly complicated form. For example, you can't have sin(w1 * t) + sin(w2 * t) or this won't work.
3. You know what the frequency and time variables are named.
4. You have a decent computer algebra system and all the standard functions that go with it.

The algorithm, under these assumptions, is:

1. Take the Laplace transform. In practice this is trivial, since most of the time when you run into problems like this your answer is originally in the Laplace domain and you end up converting back to the time domain to get your expression with the sine wave. Even if not, a decent CAS will already have a mature Laplace transform function.
2. Divide by the Laplace transform of sin(w*t). (Assuming you're using w for frequency and t for time). This gives you the transfer function, i.e. a Laplace-domain constant that an input sine wave gets multiplied by.
3. As t -> infinity, i.e. for signals that have existed for infinite time, the Laplace variable s is equivalent to i*w, so by simply substituting i*w for s, you get the transfer function in the frequency domain, in effect the Fourier transform of the transfer function.
4. For a sine wave of frequency w, the gain is the complex absolute value of the transfer function. (The complex absolute value is the distance of some complex number from the origin in the complex plane, sqrt( realpart( my_number)^2 + imagpart( my_number)^2).
5. The phase offset is just arctan( imaginarypart(transferFunction) / realpart( transferFunction)).