Sine-Tone generator produces undesired harmonics

Question

I am writing a program to generate sinewave audio tones at 44.1khz sampling rate. The 2 plots below contain 720 audio samples on the X axis and a sine wave that goes between +/-4000 on the Y axis (4000 represents the amplitude/volume in this case).

The contents of the 720 sample sine wave (1,440 bytes of PCM 16 signed audio) are sent as back to back datagrams to an audio application where it should reproduce a pure test tone. In order for a pure tone, I need to produce the lower plot (I have to adjust the frequency by hand until I get a continuous plot).

My questions is how do I mathematically calculate the set of valid frequencies to fit the above.

This is how I produced the above plots:

// check sampleRate, channelCount and sampleFormat
if (rAudioFormat.isValid()) {
    const auto sampleFormat = rAudioFormat.sampleFormat();
    const auto channelCount = rAudioFormat.channelCount();
    const auto bytesPerFrame = rAudioFormat.bytesPerFrame();
    const auto volume = aVolumePercent / 100.0f;
    // payload must have room for at least 1 frame
    if (aAudioSizeBytes >= bytesPerFrame) {
        static std::default_random_engine gRandomNumberGenerator;
        // each frame contains inputChannelCount audio samples
        const qint64 frameCount = aAudioSizeBytes / bytesPerFrame;
        // pad size in bytes at the end of the datagram.
        const qint64 paddingBytes = aAudioSizeBytes % bytesPerFrame;
        // allocate raw data for an unpadded audio sample
        const auto payload = std::make_unique<char[]>(aAudioSizeBytes - paddingBytes);
        switch (aSignalGenerator) {
        case SignalGenerator::Monotone:
            {
                if ((!rGeneratorParams.empty()) &&
                    (static_cast<size_t>(rAudioFormat.sampleRate()) > 0u)) {
                    const auto frequency = rGeneratorParams[0];
                    for (auto t = 0; t < frameCount; ++t) {
                        const auto value = qSin(
                            2 * M_PI * frequency * (t % frameCount) /
                            rAudioFormat.sampleRate());
                        switch (sampleFormat) {
                        . . .
                        case QAudioFormat::SampleFormat::Int16:
                            {
                                constexpr auto gMaxVal = std::numeric_limits<int16_t>::max();
                                const auto dest = reinterpret_cast<int16_t*>(payload.get());
                                for (auto i = 0; i < channelCount; i++) {
                                    dest[t * channelCount + i] =
                                        static_cast<int16_t>(gMaxVal * volume * value);
                                }
                            }
                            break;

So far, I have come up with the following bit of code which does part of the job (I think). I know I need get the even (as it needs to fit a sine) divisors of 720, but I don't quite know how to complete this code to return the frequency values corresponding to each of the divisors. Given the sampling rate, I should probably also consider nyquist to limit the maximum reproducible tone frequency.

    //! helper function to calculate produce pure tones
    //! given nTotalSamples with given nSampleRate.
    std::vector<int> getDivisors(int nSampleRate, int nTotalSamples, int nearestMin) {
        std::vector<int> divisors;
        if ((nSampleRate> 0) && (nearestMin > 0) &&
            (nearestMin <= nTotalSamples)) {
            for (int i = nearestMin; i < nTotalSamples; ++i) {
                if ((nTotalSamples % i) == 0) {
                    // calculate the corresponding tone frequency
                    divisors.emplace_back(i);
                }
            }
        }
        return divisors;
    }

Answer 1

Your frame duration is number of samples divided by sample rate. The number of oscillations is frequency times duration. You want that to be an integer. So simply compute the number you would get, then round that to the nearest integer and convert back to a frequency:

double frameDuration = samplesPerFrame / sampleRate;
double fractionalOscillations = frameDuration * desiredFrequency;
double wholeOscillations = round(fractionalOscillations);
double goodFrequency = wholeOscillations / frameDuration;

Note I'm using double so I don't have to worry about any rounding effects along the way. This could lead to each individual oscillation taking a fractional number of samples. The way you compute your individual values – by just feeding an argument to a sine function – I would not expect this to be a problem.

Also note that you have to use compatible units. If you use sample rate as samples per second, and frequency as Hz, that works fine since both have unit 1/s.

Answer 2

with continuous signal synthesis there are 2 major problems:

1. aliasing

   in case your f frequency is not exact integer divisor of sampling frequency fsampling the output will alias to 2 frequencies one slightly smaller and one slightly bigger than your desired frequency and they will alternate ...

   to avoid this you can round output frequency to integer divisor of sampling frequency so:

   f' = fsampling/floor(fsampling/f) // or ceil or round
   

2. fixed (finite) buffer size

   so you have a 720 samples in your buffer but you have to use it as infinite (ring FIFO) buffer so if periods of your signal is not exact divisor of buffer size you have to remember the phase of last sample and continue next buffer from that point ...

   see: Find start point (time) of each cycle in a sine wave

   If you don't then it will create glitches on joint with the next buffer producing usually many many more (even higher) frequencies...

   In case you are sending just single buffer that should loop itself then you have to round your frequency f period T (in the same manner as above) to be integer divisor of buffer duration period Tbuffer...

   Tbuffer = samples / fsampling
   T = 1/f
   T' = Tbuffer/floor(Tbuffer/T) // or ceil or round
   f' = 1/T'
   

   this will also avoid the alias problem from previous bullet as Tbuffer is multiple of Tsampling...