# Matthew van Eerde's web log

• #### Generating a sine wave with phase

I alluded to a bug in my last post but didn't explicitly say what it was.  Hint: it's all about phase.

So let's talk about phase for a minute.  Sine waves are periodic... although they are frequently plotted over very large intervals, their natural domain is a single period.

(original)

Exercise - The original version, when printed, makes a nice armband.

Advanced Exercise - generate a one-and-a-half period sine wave, print it on both sides of a piece of paper, and make a Möbius strip.

Well, that's a step, but it's not really getting at the nature of the sine wave.  Yes, we notice that the beginning value (0) is the same as the ending value... and the beginning slope (+1) is the same as the ending slope... and there is, perhaps, a certain smoothness beyond even slope (in fact, all derivatives agree at the crossover point...)

If you see a fact, try to see it as intuitively as possible.

Do the first exercise above - print out the original image, cut it out along the lines, and make a cylinder by connecting the short ends together.  Rotate the cylinder through various orientations, paying particular attention to the location of the sine wave at all times.

What you get is hard to render in a two-dimensional blog post, but here's a stab at it (Matlab script below:)

(original)

Notice that this representation of the sine wave is an ellipse.

Exercise - Prove that the red line in the above graph is planar (Hint: you can type the equation of the plane in three keystrokes.)

Advanced Exercise - Prove that it's an ellipse.

Advanced Exercise - Make a Keplerian observation. (Hint: note the angles at which the red line crosses the dotted black lines.)

So what does this have to do with phase?

Everything.

There's another hint in the documentation of double sin(double), which I also linked in my last post:

sin returns the sine of x [phase]. If x [phase] is greater than or equal to 263, or less than or equal to –263, a loss of significance in the result occurs. (my emphasis)

A final hint: again from my last post, "as a side effect, we get a kinda-sorta implementation of ϕ [phase]."

When you find yourself with a kinda-sorta implementation of something, that is occasionally a hint from above that your life will be easier if you do some self-examination... perhaps a full-blown implementation of that thing will be (a) easier and (b) more what you actually want.

No more hints.  If you want to figure it out for yourself, stop reading now.

OK, here's the scoop.  The thing you feed to sin(...) is a radian measurement.  This is naturally between 0 and 2π (or, perhaps, -π to π if you prefer.)  You can feed it things that are out of this range... say, sin(100.0)... but this is a little out of the natural and there are consequences.

double-type numbers, and floating point numbers in general, can represent very small numbers to a very high degree of precision... and they can represent very large numbers with the same relative degree of precision.  See IEEE754 for gory details.

An illustration is in order.  Richard Feynman, celebrity physicist, tells a story of his youth.  He claimed to be able to solve, in one minute, to within 10% accuracy, any problem that can be stated in ten seconds.

He did pretty well until someone came up with "tan(10100)".

If you feed bigger and bigger numbers to sin(...), eventually it starts giving you sloppy answers.  It has to.  It's not getting enough information to give you a good answer.

I was able to test the tone generation function and measure the signal-to-noise ratio of the tones it was producing after varying lengths of sample time.  It started off very well, but after very long times (days, weeks, years) the SNR would get slowly compromised until, after millions of years, the signal disappeared.

That is not the bug.

But it is a serious clue to the bug.

If the worst-case scenario for an implementation is that the signal will depreciate over a million years, ship it.

Take your cylinder that you made above.  Place it in front of you so that the taped-together part is directly in front of your nose, and the part to the right has the red line going up.  This is a sine wave.

Rotate it clockwise 90 degrees, so the part in front of your nose has the red line brushing the top.  Now it is a cosine wave.

Rotate it counter-clockwise 60 degrees, so the part in front of your nose has the red line crossing the dotted line.  Now it is a sine wave with a known phase (30 degrees.)

In fact, rotate it any angle (say, ϕ) and it's a sine wave with known phase ϕ.

We can't do that with our tone generator.

That is the bug.  We're implementing a tweaked version of phase which has a forgivable side effect, but which is not a full implementation of phase.

Why do we care about implementing phase?

Well, for example, one problem I sometimes have when testing audio fidelity is I connect the stereo cables to the Audio Precision equipment the wrong way (left-to-right and right-to-left.)  If I could generate the test tones out of phase by 90 degrees, with the left channel in front of the right channel, I could take a phase measurement and log a "you have the cables backwards" error.  I send different frequencies over the cables, so the ifFrameOffset parameter is of little use here.

So what's the improved implementation?

// given a buffer, fills it with sine wave data
HRESULT FillBufferWithSine(
LPCWAVEFORMATEX pWfx, // includes sample rate
double lfAmplitude,
double lfFrequency,
UINT64 ifFrameOffset,
PBYTE pBuffer,
UINT32 cbBytes
);

do

// given a buffer, fills it with sine wave data
HRESULT FillBufferWithSine(
LPCWAVEFORMATEX pWfx, // includes sample rate
double lfAmplitude,
double lfFrequency,

// in/out
//     On input, initial phase (0 to 2pi).
//     On output, phase of next sample.
double *plfPhase,

PBYTE pBuffer, // out - byte count is cbBytes
UINT32 cbBytes
);

Now, if you want to generate a lot of contiguous sine tone buffers, you do something like this:

double lfPhase = 0.0; // start at 0... or somewhere else
while (...) {
// ...
// FillBufferWithSine takes care of updating lfPhase as necessary
hr = FillBufferWithSine(pWfx, lfAmp, lfFreq, &lfPhase, pBuffer, cbBytes);
// ...
}

As a side effect of this new phase feature, the quality degradation non-issue goes away.  So the generated tones will sound good forever... which is kind of neat... but not the motivation for this fix.

Why did this happen?

Because there are many ways to measure time.  You can count samples... you can count milliseconds... you can count hundred-nanosecond intervals... you can count CPU flops... you can do so many different things.  Some are more natural than others in different contexts.

The reason this happened is because a time measure that was very natural in one context (counting frames: playing an audio stream) was somewhat artificially shoved into another context (generating a sine wave) where there was another, more natural context (changing phase.)  Alas, context switching is a frequent source of bugs.

Matlab script to generate the three-dimensional rendering of a sine wave:

t = 0:pi/50:2*pi;
h = -1:0.5:1;
z = ones(size(t))' * h;

sinx = cos(t);
siny = sin(t);
sinz = sin(t);

plot3( ...
sinx, siny, z, 'k:', ...
sinx, siny, sinz, 'r' ...
);

• #### Generating a sine wave

One common task on the Windows Sound test team is taking a known-good sine wave, passing it through some chunk of code, and looking at what comes out.  Is it a sine wave?  Is it the same frequency that it should be?  Etc., etc.

So the quality of the sine wave generator is important.

There's a commonly available sine wave generator: double sin(double).   It's fairly simple to use this to generate a buffer, if you keep a few things in mind.  But there are pitfalls too.

In full generality, a sine wave can be completely described by four measurements:

1. Amplitude: half the difference between the maximum and minimum samples.
There's a mathematical reason to use half the difference instead of the whole difference.  We'll go into this later.
Yes, there were some bugs due to this "use half the difference" convention.
For our purposes, we'll consider "full scale" to be from -1.0 to 1.0; the amplitude of a full-scale sine wave is 1.0.  the amplitude of a -3 dB FS sine wave is 10-3/20, or about 0.707945...
2. Frequency or Period.  Frequency is the number of complete cycles in a fixed unit of time (usually one second.)  Period is the time length of a single complete cycle.  Frequency * Period = 1, unless one of them is zero (in which case the other is infinite.)
3. dc.  dc is the average value over a complete cycle.
4. Phase.  This tells where in the cycle the first sample is.  If you measure the angle in radians... which is usual... then the phase is a number between 0 and 2π.

These can all be linked together by a single formula which tells you the mathematical value of a sine wave with these features at any given time.  Let:

a = amplitude (0.0 to 1.0)
f = frequency (in cycles/sec)
c = dc (0.0 to 1.0; to avoid clipping, 0.0 to (1.0 - a))
ϕ = starting phase (0.0 to 2π)
t = the time of interest in seconds

Then

w(t) = a sin(2πft + ϕ) + c

Note in passing that w(0) = a sin((2πf)0 + ϕ) + c = a sin(ϕ) + c is independent of f.

So much for mathematics.  When you want to realize a sine wave in a PCM wave format, you have to discretize the time dimension and the sample dimension.  To discretize the time dimension you need a sample rate (usually 44100 samples per second) and to discretize the sample dimension you need to specify a sample format (usually 16-bit int).

Let's talk about discretizing the time dimension first.  Let s = the sample rate, and i be the sample of interest.  Now:

w(i) = a sin(2πfi / s + ϕ) + c

As a design simplification, we collect constants, and eliminate ϕ and c...

w(i) = a sin((2πf / s) i)

And we have an implementable function

// given a buffer, fills it with sine wave data
HRESULT FillBufferWithSine(
LPCWAVEFORMATEX pWfx, // includes sample rate
double lfAmplitude,
double lfFrequency,
PBYTE pBuffer,
UINT32 cbBytes
);

This would start at i = 0, loop through i = cbBytes / pWfx->nBlockAlign, and set sample values for each i.

It is often necessary to get continuous sine wave data in chunks... say, 1 KB at a time.  It is important that the phase of the sine wave match at chunk boundaries.  This suggests another parameter:

// given a buffer, fills it with sine wave data
HRESULT FillBufferWithSine(
LPCWAVEFORMATEX pWfx, // includes sample rate
double lfAmplitude,
double lfFrequency,
UINT32 ifFrameOffset,
PBYTE pBuffer,
UINT32 cbBytes
);

... where ifFrameOffset is the count of frames to skip.  Now a client can fill a 1 KB buffer with 256 frames of 16-bit int stereo sine wave, pass it off to someone else, fill another 1 KB buffer starting at frame offset 256 with 16-bit int stereo sine wave, etc., and the effect is an almost infinite continuous sine wave.  (Oh, and as a side effect, we get a kinda-sorta implementation of ϕ.)

Almost.

Why almost?

Well, because I used UINT32 as the frame offset.  Suppose you, as a client, use a WAVEFORMATEX with a nice high sample rate like 384 kHz.  A UINT32 wraps at 232.   At 384,000 samples per second, you'll burn through 232 samples in (232 / 384000) seconds... that's, um, 1184.8106666... seconds... three hours, six minutes, twenty-five seconds.  This is uncomfortably close to reasonable.

At that point, one of two things happen.  If you were clever in your chosen frequency, you will wrap back to a point in the sine wave where the phase is the same as what you want anyway... and there will be no observable artifact in the generated sine wave.  Or, if you were unlucky, you will wrap back to a point in the sine wave where the phase is different - and there will be a discontinuity.  This will probably result in an audible "pop" artifact in the signal tone.

Our test tone generator had this problem.  We knew about it - we didn't care.

But sometimes little things nag at me and I have to fix them or I can't sleep at night.  It's a fixation... a terrible disease... but there it is.

So I fixed it - almost.  Now we use a UINT64 for the sample offset.  That won't wrap until 264 samples.  Even at 384 kHz, that won't wrap until (264/384000) seconds - that's, um,  48038396025285.29 seconds... 1,521,800-odd years.  I'll let a tester from that era file a bug then. :-)

Almost.

Why almost?

Well, it turns out that fixing that bug uncovered another bug.  We know about it - we don't care.  (But sometimes little things nag at me...)

Exercise... what is the bug?

Page 1 of 1 (2 items)