Aliasing of oscillator

pulse.cpp

// Square Wave Oscillator Comparison
// This single-file C++ program generates "out_<method>_<sampleRate>.wav".
//   M8P2
//     Original oscillator of Magical 8-bit Plug 2
//   Sines
//     Sum of sines below sampleRate / 2
//   ManySines
//     Sum of sines below 2.5 MHz (including sines above sampleRate / 2)

#define _USE_MATH_DEFINES
#include <stdio.h>
#include <cmath>
#include <vector>
#include <string>
#include <sstream>

const size_t kSampleRate = 44100;             // Sample rate
const size_t kSampleLength = kSampleRate * 2; // 2 sec
const double kMinNoteNumber = 48;             // C3
const double kMaxNoteNumber = 132;            // C10

///// Oscillators

class AbstractOscillator
{
public:
    virtual std::string getName() const { return std::string("Abstract"); }
    virtual double getNextSample(double frequency) { return 0.0; }
    virtual ~AbstractOscillator() {}

protected:
    AbstractOscillator() {}

private:
    AbstractOscillator(const AbstractOscillator &);
    AbstractOscillator &operator=(const AbstractOscillator &);
};

class OscillatorM8P2 : public AbstractOscillator
{
private:
    double _currentAngle = 0.0;

public:
    std::string getName() const { return std::string("M8P2"); }

    double getNextSample(double frequency)
    {
        double rate = 1.0; // duty = 50%
        double output = _currentAngle < rate * M_PI ? -1.0 : 1.0;

        _currentAngle += 2.0 * M_PI * frequency / kSampleRate;
        while (_currentAngle > 2.0 * M_PI)
            _currentAngle -= 2.0 * M_PI;

        return output;
    }
};

class OscillatorSines : public AbstractOscillator
{
private:
    double _current = 0.0;

public:
    std::string getName() const { return std::string("Sines"); }

    double getNextSample(double frequency)
    {
        // add sines below the kSampleRate/2
        double output = 0.0;
        double currentAngle = 2.0 * M_PI * _current;
        double c = 1.0;
        while (c * frequency <= kSampleRate / 2.0)
        {
            output -= std::sin(currentAngle * c) / c;
            c += 2.0;
        }
        output *= 4.0 / M_PI;

        _current += frequency / kSampleRate;
        _current -= std::floor(_current);
        return output;
    }
};

class OscillatorManySines : public AbstractOscillator
{
private:
    double _current = 0.0;

public:
    std::string getName() const { return std::string("ManySines"); }

    double getNextSample(double frequency)
    {
        const double f = 2500000.0;

        // add sines below the f
        double output = 0.0;
        double currentAngle = 2.0 * M_PI * _current;
        double c = 1.0;
        while (c * frequency <= f)
        {
            output -= std::sin(currentAngle * c) / c;
            c += 2.0;
        }
        output *= 4.0 / M_PI;

        _current += frequency / kSampleRate;
        _current -= std::floor(_current);
        return output;
    }
};

///// Utilities

// frequency sweep
double getSweepFrequency(double current, double maximum)
{
    double t = current / maximum;
    double note = t * (kMaxNoteNumber - kMinNoteNumber) + kMinNoteNumber;
    return 440.0 * std::exp2((note - 69.0) / 12.0);
}

// writes WAV file
// monoral, 32-bit floating point
// (only works on little-endian machines)
void writeWaveFile(const std::string &filename, const float *waveform, size_t length, int32_t samplerate)
{
    FILE *fout = fopen(filename.c_str(), "wb");
    if (fout == NULL)
    {
        printf("  Error: Unable to write WAV file: '%s'\n", filename.c_str());
        return;
    }

    uint8_t header[44] = {};
    uint16_t channels = 1;
    uint16_t samplesize = sizeof(*waveform);
    uint16_t blocksize = samplesize * channels;
    uint32_t filesize = sizeof(header) + length * blocksize;
    *((uint32_t *)header) = 0x46464952;
    *((uint32_t *)(header + 4)) = filesize - 8;
    *((uint32_t *)(header + 8)) = 0x45564157;
    *((uint32_t *)(header + 12)) = 0x20746D66;
    *((uint32_t *)(header + 16)) = 16;
    *((uint16_t *)(header + 20)) = 3;
    *((uint16_t *)(header + 22)) = channels;
    *((uint32_t *)(header + 24)) = samplerate;
    *((uint32_t *)(header + 28)) = samplerate * blocksize;
    *((uint16_t *)(header + 32)) = blocksize;
    *((uint16_t *)(header + 34)) = samplesize * 8;
    *((uint32_t *)(header + 36)) = 0x61746164;
    *((uint32_t *)(header + 40)) = filesize - sizeof(header);
    fwrite(header, 1, sizeof(header), fout);
    fwrite(waveform, samplesize, length, fout);
    fclose(fout);

    printf("  Wrote '%s'\n", filename.c_str());
}

///// Main

int main()
{
    // list of the oscillator instances
    std::vector<AbstractOscillator *> oscillators = {
        new OscillatorM8P2(),
        new OscillatorSines(),
        new OscillatorManySines(),
    };

    puts("Settings:");
    printf("  Sample rate:       %9.2lf Hz\n", (double)kSampleRate);
    printf("  Sample rate / 2:   %9.2lf Hz\n", (double)kSampleRate / 2.0);
    printf("  Minimum frequency: %9.2lf Hz\n", getSweepFrequency(0, kSampleLength - 1));
    printf("  Maximum frequency: %9.2lf Hz\n", getSweepFrequency(kSampleLength - 1, kSampleLength - 1));

    printf("\nGenerating square-wave with %d types of oscillators...\n", (int)oscillators.size());
    for (auto current = oscillators.begin(); current != oscillators.end(); current++)
    {
        auto osc = *current;

        float waveform[kSampleLength];
        for (size_t i = 0; i < kSampleLength; i++)
        {
            double frequency = getSweepFrequency(i, kSampleLength - 1);
            double output = 0.5 * osc->getNextSample(frequency);
            waveform[i] = (float)output;
        }

        std::stringstream ss;
        ss << "out_";
        ss << osc->getName();
        ss << "_";
        ss << kSampleRate;
        ss << ".wav";

        std::string filename = ss.str();
        writeWaveFile(filename, waveform, kSampleLength, kSampleRate);

        delete osc;
        *current = nullptr;
    }

    puts("\nAll processes have done!");
    return 0;
}