Loss of double precision C#

Question

I am experiencing precision loss when im using doubles and i cant seem to find where the precision is lost. I am writing a software synthesizer and for some reason the input frequency for an oscillator(sound-wave generator) gets heavily aliased somewhere along the way. The oscillator is supposed to generate a triangle waveform and the correct frequency is passed to the method but when the result is played in the speakers i can clearly hear that the sound snaps to different frequencies.

I am using NAudio (http://naudio.codeplex.com/) and this method is run once for every sample i want to generate.

Here is the code:

double counter = 0;
int samplecounter = 0;

public double GetNextTriangle(double frequency)
{
    double samplesperwave = (double)Parent.WaveFormat.SampleRate / frequency;
    double length = (double)samplecounter / samplesperwave;
    if (length.CompareTo(0.25) == -1)
    {
        counter = length * 4.0;
    }
    else if (length.CompareTo(0.75) == -1)
    {
        counter = 1.0 - (length - 0.25) * 4.0;
    }
    else
    {
        counter = -1.0 + (length - 0.75) * 4.0;
    }
    samplecounter++;
    samplecounter = samplecounter > samplesperwave ? 0 : samplecounter;
    return counter;
}

Thanks in advance! //Mats

Answer 1

Your problem is not one of precision. The issue is that your function does not define a triangle wave. Each time samplecounter is reset to 0, the next returned value from the function is 0. Next time round length is set to 0, the first if branch is executed and counter is set to 0.

When faced with a problem like this you should plot the output. Had you done so you would have seen immediately that your function does not produce a triangle wave form. Here's some code that does:

static double GetNextTriangle(double frequency)
{
    double samplesperwave = Parent.WaveFormat.SampleRate / frequency;
    double t = samplecounter / samplesperwave;
    counter = 1.0 - 4.0 * Math.Abs(Math.Round(t - 0.25) - (t - 0.25));
    samplecounter++;
    return counter;
}

Again, I cannot stress enough that you must plot and visualize the output from your code to gain insight into its behaviour.

Answer 2

To try and get a better idea of what's wrong, let's visualize the data. Here's two cycles of a 2 Hz wave at 41 samples per second:

There's a blip around sample 20. It looks like the samples aren't fitting a triangle wave properly.

Let's go over the math of a triangle wave. What you're doing is sampling a continuous function at discrete points. First, define a continuous triangle wave of frequency f (or 1/T) in hertz:

tri(t) =   {  4*f*t,            0     <= t < T/4   }
       =   { -4*f*(t - T/2),    T/4   <= t < 3*T/4 }
       =   {  4*f*(t - T),      3*T/4 <= t < T     }
       =   tri(t - T)  [it's periodic!]

You now want to sample this continuous function. So you define a sampling rate, s (or 1/U) in samples per second. Now, the nth sample will simply be tri(n*U):

tri[n] = tri(n*U)
       =   {  4*f*n*U,            0     <= n*U < T/4   }
       =   { -4*f*(n*U - T/2),    T/4   <= n*U < 3*T/4 }
       =   {  4*f*(n*U - T),      3*T/4 <= n*U < T     }

Let's clean it up a bit by defining the normalized period, P = T/U, and the normalized frequency, F = f/s = U/T:

tri[n] =   {  4*F*n,            0     <= n < P/4   }
       =   { -4*F*(n - P/2),    P/4   <= n < 3*P/4 }
       =   {  4*F*(n - P),      3*P/4 <= n < P     }

Now we get this:

Don't worry about the apparent "blips" at the tips; they are to be expected and you shouldn't try to avoid them.

Here's the code:

public static double GetNextTriangle(int sample, double frequency, double sampleRate)
{
    double T = 1d / frequency;
    double U = 1d / sampleRate;

    double P = T / U;
    double F = U / T;

    double n = (double)sample;
    n %= P; // restrict n to the domain [0, P)

    if ((n >= 0) && (n < (P / 4d)))
    {
        return 4d * F * n;
    }
    else if ((n >= (P / 4d)) && (n < (3d * P / 4d)))
    {
        return -4d * F * (n - (P / 2d));
    }
    else // if ((n >= (3d * P / 4d)) && (n < P))
    {
        return 4d * F * (n - P);
    }
}

Answer 3

The floating point operations here are well conditionned. But you may see an alisaing problem related to sampling if 1/frequency is not a multiple of 1/sampleRate.

Try this matlab code if you can

sampleRate=5000;
frequency=700;
sampleperwave=sampleRate/frequency;
samplecounter=0:floor(sampleperwave);
samplecounter=repmat(samplecounter,1,5);
length=samplecounter/sampleperwave;
wave=-1+4*(length-0.75);
wave(length<0.75)=1-4*(length(length<0.75)-0.25);
wave(length<0.25)=4*length(length<0.25);
figure; stem(wave); hold on; plot(wave,'r')

You could try to declare samplecounter as double and increment it with a modulo

samplecounter++;
samplecounter = samplecounter % samplesperwave;

Or back in matlab

samplecounter=0:length(samplecounter)-1;
samplecounter=rem(samplecounter,sampleperwave);
len=samplecounter/sampleperwave;
wave=-1+4*(len-0.75);
wave(len<0.75)=1-4*(len(len<0.75)-0.25);
wave(len<0.25)=4*len(len<0.25);
figure; stem(wave); hold on; plot(wave,'r')

Answer 4

The whole approach seems wrong. I'd go with something like this:

public class RectWave
{ // triangular wave in range [0,1] with given frequency
    public double GetNextSample()
    {
        if (_ptr == _wave.Length)
            _ptr = 0;
        return _wave[_ptr++];
    }

    public RectWave(int sampleRate = 441000, int period = 20)
    { 
        _sampleRate = sampleRate;
        _period = period;
        _wave = new double[(int)(_sampleRate * _period / 1000.0)];  // one cycle
        _ptr = 0;
        BuildWave();
    }

    private void BuildWave()
    {
        double dt = 1000.0 / _sampleRate;   // in msec
        double slope = 1.0 / (_period / 2.0);   // 1.0 => peek 
        for (int i=0; i <= _wave.Length/2; ++i)
        {
            _wave[i] = _wave[_wave.Length - i - 1] = i * dt * slope; 
        }
    }

    private int _sampleRate;    // in Hz
    private int _period;        // in msec
    private double[] _wave;     // the data
    private int _ptr;           // current sample
}