`index2freq` Simplification

[ryukau]

Continued from: https://github.com/ryukau/LV2Plugins/issues/29

The idea is to do computation on normalized frequency.

Current index2freq implementation.

index2freq(index) = ((index-index')*ma.SR) : ba.sAndH(abs(index-index')<0.5);

The output of index2freq is (index - index') * ma.SR with sample and hold (S&H). I'll omit S&H section later, because it doesn't change the result of addition/multiplication (only if there's no delay involved, and that's the case here).

indexAA

In indexAA, there's a following expression:

(index2freq(fund)-(ma.SR/256)) / (ma.SR/8)

Expanding index2freq:

((fund - fund') * ma.SR - ma.SR / 256) / (ma.SR/8)

This can be simplified to:

(8/ma.SR) * (fund - fund') * ma.SR - (8/ma.SR) * ma.SR / 256

8 * (fund - fund') - 8 / 256

resAA

resAA expression:

res * index2freq(fund) : max(0) : min(ma.SR/4) / index2freq(fund);

If I understand : operator correctly, this can be rewritten to:

min(
  max(
    res * index2freq(fund), // a
    0                       // b
  ),
  ma.SR/4                   // c
) / index2freq(fund);

Applying division of index2freq(fund) to line a, b, c:

min(
  max(
    res, // a
    0    // b
  ),
  (ma.SR/4) / index2freq(fund) // c
);

Expanding index2freq:

min(
  max(res, 0),
  (ma.SR/4) / ((fund - fund') * ma.SR)
);

ma.SR can be eliminated:

min(
  max(res, 0),
  (1/4) / (fund - fund')
);

Rewriting this with : operator:

res : max(0) : min(1 / (4 * (fund - fund')));

[magnetophon]

I suspected faust would do this simplification for me, so I tested it, and it indeed does, up to a certain point.

Without the sAndH, the cpp is almost identical:

This faust code:

import("stdfaust.lib");
process =
(index2freq(fund)-(ma.SR/256)) / (ma.SR/8);
fund = os.lf_sawpos(440);
index2freq(index)        = ((index-index')*ma.SR);

Gives me this cpp:

/* ------------------------------------------------------------
name: "test"
Code generated with Faust 2.27.1 (https://faust.grame.fr)
Compilation options: -lang cpp -scal -ftz 0
------------------------------------------------------------ */

#ifndef  __mydsp_H__
#define  __mydsp_H__

#ifndef FAUSTFLOAT
#define FAUSTFLOAT float
#endif 

#include <algorithm>
#include <cmath>
#include <math.h>

#ifndef FAUSTCLASS 
#define FAUSTCLASS mydsp
#endif

#ifdef __APPLE__ 
#define exp10f __exp10f
#define exp10 __exp10
#endif

class mydsp : public dsp {

 private:

    int fSampleRate;
    float fConst0;
    float fRec0[2];

 public:

    void metadata(Meta* m) { 
        m->declare("filename", "test.dsp");
        m->declare("maths.lib/author", "GRAME");
        m->declare("maths.lib/copyright", "GRAME");
        m->declare("maths.lib/license", "LGPL with exception");
        m->declare("maths.lib/name", "Faust Math Library");
        m->declare("maths.lib/version", "2.3");
        m->declare("name", "test");
        m->declare("oscillators.lib/name", "Faust Oscillator Library");
        m->declare("oscillators.lib/version", "0.1");
        m->declare("platform.lib/name", "Generic Platform Library");
        m->declare("platform.lib/version", "0.1");
    }

    virtual int getNumInputs() {
        return 0;
    }
    virtual int getNumOutputs() {
        return 1;
    }
    virtual int getInputRate(int channel) {
        int rate;
        switch ((channel)) {
            default: {
                rate = -1;
                break;
            }
        }
        return rate;
    }
    virtual int getOutputRate(int channel) {
        int rate;
        switch ((channel)) {
            case 0: {
                rate = 1;
                break;
            }
            default: {
                rate = -1;
                break;
            }
        }
        return rate;
    }

    static void classInit(int sample_rate) {
    }

    virtual void instanceConstants(int sample_rate) {
        fSampleRate = sample_rate;
        fConst0 = (440.0f / std::min<float>(192000.0f, std::max<float>(1.0f, float(fSampleRate))));
    }

    virtual void instanceResetUserInterface() {
    }

    virtual void instanceClear() {
        for (int l0 = 0; (l0 < 2); l0 = (l0 + 1)) {
            fRec0[l0] = 0.0f;
        }
    }

    virtual void init(int sample_rate) {
        classInit(sample_rate);
        instanceInit(sample_rate);
    }
    virtual void instanceInit(int sample_rate) {
        instanceConstants(sample_rate);
        instanceResetUserInterface();
        instanceClear();
    }

    virtual mydsp* clone() {
        return new mydsp();
    }

    virtual int getSampleRate() {
        return fSampleRate;
    }

    virtual void buildUserInterface(UI* ui_interface) {
        ui_interface->openVerticalBox("test");
        ui_interface->closeBox();
    }

    virtual void compute(int count, FAUSTFLOAT** inputs, FAUSTFLOAT** outputs) {
        FAUSTFLOAT* output0 = outputs[0];
        for (int i = 0; (i < count); i = (i + 1)) {
            fRec0[0] = (fConst0 + (fRec0[1] - std::floor((fConst0 + fRec0[1]))));
            output0[i] = FAUSTFLOAT((8.0f * ((fRec0[0] - fRec0[1]) + -0.00390625f)));
            fRec0[1] = fRec0[0];
        }
    }

};

#endif

And this faust code:

import("stdfaust.lib");
process =
  8 * (fund - fund') - 8 / 256;

Gives me this diff in cpp:

<           output0[i] = FAUSTFLOAT((8.0f * ((fRec0[0] - fRec0[1]) + -0.00390625f)));
---
>           output0[i] = FAUSTFLOAT(((8.0f * (fRec0[0] - fRec0[1])) + -0.03125f));

With the sAndH, I get different results:

This faust code:

import("stdfaust.lib");
process =
  (index2freq(fund)-(ma.SR/256)) / (ma.SR/8);
fund = os.lf_sawpos(440);
index2freq(index)        = ((index-index')*ma.SR) : ba.sAndH(abs(index-index')<0.5);

gives me this cpp:

/* ------------------------------------------------------------
name: "test"
Code generated with Faust 2.27.1 (https://faust.grame.fr)
Compilation options: -lang cpp -scal -ftz 0
------------------------------------------------------------ */

#ifndef  __mydsp_H__
#define  __mydsp_H__

#ifndef FAUSTFLOAT
#define FAUSTFLOAT float
#endif 

#include <algorithm>
#include <cmath>
#include <math.h>

#ifndef FAUSTCLASS 
#define FAUSTCLASS mydsp
#endif

#ifdef __APPLE__ 
#define exp10f __exp10f
#define exp10 __exp10
#endif

class mydsp : public dsp {

 private:

    int fSampleRate;
    float fConst0;
    float fConst1;
    float fConst2;
    float fRec1[2];
    float fRec0[2];
    float fConst3;

 public:

    void metadata(Meta* m) { 
        m->declare("basics.lib/name", "Faust Basic Element Library");
        m->declare("basics.lib/version", "0.1");
        m->declare("filename", "test.dsp");
        m->declare("maths.lib/author", "GRAME");
        m->declare("maths.lib/copyright", "GRAME");
        m->declare("maths.lib/license", "LGPL with exception");
        m->declare("maths.lib/name", "Faust Math Library");
        m->declare("maths.lib/version", "2.3");
        m->declare("name", "test");
        m->declare("oscillators.lib/name", "Faust Oscillator Library");
        m->declare("oscillators.lib/version", "0.1");
        m->declare("platform.lib/name", "Generic Platform Library");
        m->declare("platform.lib/version", "0.1");
    }

    virtual int getNumInputs() {
        return 0;
    }
    virtual int getNumOutputs() {
        return 1;
    }
    virtual int getInputRate(int channel) {
        int rate;
        switch ((channel)) {
            default: {
                rate = -1;
                break;
            }
        }
        return rate;
    }
    virtual int getOutputRate(int channel) {
        int rate;
        switch ((channel)) {
            case 0: {
                rate = 1;
                break;
            }
            default: {
                rate = -1;
                break;
            }
        }
        return rate;
    }

    static void classInit(int sample_rate) {
    }

    virtual void instanceConstants(int sample_rate) {
        fSampleRate = sample_rate;
        fConst0 = std::min<float>(192000.0f, std::max<float>(1.0f, float(fSampleRate)));
        fConst1 = (8.0f / fConst0);
        fConst2 = (440.0f / fConst0);
        fConst3 = (0.00390625f * fConst0);
    }

    virtual void instanceResetUserInterface() {
    }

    virtual void instanceClear() {
        for (int l0 = 0; (l0 < 2); l0 = (l0 + 1)) {
            fRec1[l0] = 0.0f;
        }
        for (int l1 = 0; (l1 < 2); l1 = (l1 + 1)) {
            fRec0[l1] = 0.0f;
        }
    }

    virtual void init(int sample_rate) {
        classInit(sample_rate);
        instanceInit(sample_rate);
    }
    virtual void instanceInit(int sample_rate) {
        instanceConstants(sample_rate);
        instanceResetUserInterface();
        instanceClear();
    }

    virtual mydsp* clone() {
        return new mydsp();
    }

    virtual int getSampleRate() {
        return fSampleRate;
    }

    virtual void buildUserInterface(UI* ui_interface) {
        ui_interface->openVerticalBox("test");
        ui_interface->closeBox();
    }

    virtual void compute(int count, FAUSTFLOAT** inputs, FAUSTFLOAT** outputs) {
        FAUSTFLOAT* output0 = outputs[0];
        for (int i = 0; (i < count); i = (i + 1)) {
            fRec1[0] = (fConst2 + (fRec1[1] - std::floor((fConst2 + fRec1[1]))));
            float fTemp0 = (fRec1[0] - fRec1[1]);
            fRec0[0] = ((std::fabs(fTemp0) < 0.5f) ? (fConst0 * fTemp0) : fRec0[1]);
            output0[i] = FAUSTFLOAT((fConst1 * (fRec0[0] - fConst3)));
            fRec1[1] = fRec1[0];
            fRec0[1] = fRec0[0];
        }
    }

};

#endif

and this faust code:

import("stdfaust.lib");
process =
  8 * (fund - fund'): ba.sAndH(abs(fund-fund')<0.5) - 8 / 256;
fund = os.lf_sawpos(440);
index2freq(index)        = ((index-index')*ma.SR) : ba.sAndH(abs(index-index')<0.5);

gives me this cpp:

/* ------------------------------------------------------------
name: "test"
Code generated with Faust 2.27.1 (https://faust.grame.fr)
Compilation options: -lang cpp -scal -ftz 0
------------------------------------------------------------ */

#ifndef  __mydsp_H__
#define  __mydsp_H__

#ifndef FAUSTFLOAT
#define FAUSTFLOAT float
#endif 

#include <algorithm>
#include <cmath>
#include <math.h>

#ifndef FAUSTCLASS 
#define FAUSTCLASS mydsp
#endif

#ifdef __APPLE__ 
#define exp10f __exp10f
#define exp10 __exp10
#endif

class mydsp : public dsp {

 private:

    int fSampleRate;
    float fConst0;
    float fRec1[2];
    float fRec0[2];

 public:

    void metadata(Meta* m) { 
        m->declare("basics.lib/name", "Faust Basic Element Library");
        m->declare("basics.lib/version", "0.1");
        m->declare("filename", "test.dsp");
        m->declare("maths.lib/author", "GRAME");
        m->declare("maths.lib/copyright", "GRAME");
        m->declare("maths.lib/license", "LGPL with exception");
        m->declare("maths.lib/name", "Faust Math Library");
        m->declare("maths.lib/version", "2.3");
        m->declare("name", "test");
        m->declare("oscillators.lib/name", "Faust Oscillator Library");
        m->declare("oscillators.lib/version", "0.1");
        m->declare("platform.lib/name", "Generic Platform Library");
        m->declare("platform.lib/version", "0.1");
    }

    virtual int getNumInputs() {
        return 0;
    }
    virtual int getNumOutputs() {
        return 1;
    }
    virtual int getInputRate(int channel) {
        int rate;
        switch ((channel)) {
            default: {
                rate = -1;
                break;
            }
        }
        return rate;
    }
    virtual int getOutputRate(int channel) {
        int rate;
        switch ((channel)) {
            case 0: {
                rate = 1;
                break;
            }
            default: {
                rate = -1;
                break;
            }
        }
        return rate;
    }

    static void classInit(int sample_rate) {
    }

    virtual void instanceConstants(int sample_rate) {
        fSampleRate = sample_rate;
        fConst0 = (440.0f / std::min<float>(192000.0f, std::max<float>(1.0f, float(fSampleRate))));
    }

    virtual void instanceResetUserInterface() {
    }

    virtual void instanceClear() {
        for (int l0 = 0; (l0 < 2); l0 = (l0 + 1)) {
            fRec1[l0] = 0.0f;
        }
        for (int l1 = 0; (l1 < 2); l1 = (l1 + 1)) {
            fRec0[l1] = 0.0f;
        }
    }

    virtual void init(int sample_rate) {
        classInit(sample_rate);
        instanceInit(sample_rate);
    }
    virtual void instanceInit(int sample_rate) {
        instanceConstants(sample_rate);
        instanceResetUserInterface();
        instanceClear();
    }

    virtual mydsp* clone() {
        return new mydsp();
    }

    virtual int getSampleRate() {
        return fSampleRate;
    }

    virtual void buildUserInterface(UI* ui_interface) {
        ui_interface->openVerticalBox("test");
        ui_interface->closeBox();
    }

    virtual void compute(int count, FAUSTFLOAT** inputs, FAUSTFLOAT** outputs) {
        FAUSTFLOAT* output0 = outputs[0];
        for (int i = 0; (i < count); i = (i + 1)) {
            fRec1[0] = (fConst0 + (fRec1[1] - std::floor((fConst0 + fRec1[1]))));
            float fTemp0 = (fRec1[0] - fRec1[1]);
            fRec0[0] = ((std::fabs(fTemp0) < 0.5f) ? (8.0f * fTemp0) : fRec0[1]);
            output0[i] = FAUSTFLOAT((fRec0[0] + -0.03125f));
            fRec1[1] = fRec1[0];
            fRec0[1] = fRec0[0];
        }
    }

};

#endif

Does the difference matter much?

I didn't test res yet.

[ryukau]

The point of simplification is to make wrap around handling easier.

On sAndH version, simplified version seems faster, because it replaces fConst* variables to constant. But the difference is most likely negligible, at least if the function is only called once per frame.

I'll take a benchmark later. Currently learning faustpp and Faust benchmark tools.

[magnetophon]

What do you mean by wrap around handling?

Currently learning faustpp and Faust benchmark tools.

Cool! I want to learn those too, but I haven't gotten around to it yet. :(

[ryukau]

What do you mean by wrap around handling?

In current index2freq, index - index' can be negative if the input is normalized phase in [0, 1). However, it seems like indexAA and resAA are assuming index2freq output is always positive.

So changing from

index - index'

to

if(index < index', 1, 0) + index - index'

might improve antialiasing.

[magnetophon]

That is what : ba.sAndH(abs(index-index')<0.5) is for. IOW: when the phase wraps around, use the previous frequency.

Won't your solution have a spike in the frequency, every time the phase wraps?

[ryukau]

Ah, I missed that meaning. Sorry for that.

I now see that my code cause spike when index goes backward (decreasing).

I'll close this issue. About benchmark, I'll open new issue when I get result.