Initial commit of code.

1 files changed, 230 insertions, 0 deletions

diff --git a/src/com/jsyn/engine/MultiTable.java b/src/com/jsyn/engine/MultiTable.java
new file mode 100644
index 0000000..48b03cd
--- /dev/null
+++ b/src/com/jsyn/engine/MultiTable.java
@@ -0,0 +1,230 @@
+/*
+ * Copyright 2009 Phil Burk, Mobileer Inc
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * 
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ * 
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package com.jsyn.engine;
+
+/*
+ * Multiple tables of sawtooth data.
+ * organized by octaves below the Nyquist Rate.
+ * used to generate band-limited Sawtooth, Impulse, Pulse, Square and Triangle BL waveforms
+ * 
+ <pre>
+ Analysis of octave requirements for tables.
+
+ OctavesIndex    Frequency     Partials
+ 0               N/2  11025      1
+ 1               N/4   5512      2
+ 2               N/8   2756      4
+ 3               N/16  1378      8
+ 4               N/32   689      16
+ 5               N/64   344      32
+ 6               N/128  172      64
+ 7               N/256   86      128
+ </pre>
+ *
+ * @author Phil Burk (C) 2009 Mobileer Inc
+ */
+public class MultiTable {
+
+    public final static int NUM_TABLES = 8;
+    public final static int CYCLE_SIZE = (1 << 10);
+
+    private static MultiTable instance = new MultiTable(NUM_TABLES, CYCLE_SIZE);
+    private double phaseScalar;
+    private float[][] tables; // array of array of tables
+
+    /**************************************************************************
+     * Initialize sawtooth wavetables. Table[0] should contain a pure sine wave. Succeeding tables
+     * should have increasing numbers of partials.
+     */
+    public MultiTable(int numTables, int cycleSize) {
+        int tableSize = cycleSize + 1;
+
+        // Allocate array of arrays.
+        tables = new float[numTables][tableSize];
+
+        float[] sineTable = tables[0];
+
+        phaseScalar = (float) (cycleSize * 0.5);
+
+        /* Fill initial sine table with values for -PI to PI. */
+        for (int j = 0; j < tableSize; j++) {
+            sineTable[j] = (float) Math.sin(((((double) j) / (double) cycleSize) * Math.PI * 2.0)
+                    - Math.PI);
+        }
+
+        /*
+         * Build each table from scratch and scale partials by raised cosine* to eliminate Gibbs
+         * effect.
+         */
+        for (int i = 1; i < numTables; i++) {
+            int numPartials;
+            double kGibbs;
+            float[] table = tables[i];
+
+            /* Add together partials for this table. */
+            numPartials = 1 << i;
+            kGibbs = Math.PI / (2 * numPartials);
+            for (int k = 0; k < numPartials; k++) {
+                double ampl, cGibbs;
+                int sineIndex = 0;
+                int partial = k + 1;
+                cGibbs = Math.cos(k * kGibbs);
+                /* Calculate amplitude for Nth partial */
+                ampl = cGibbs * cGibbs / partial;
+
+                for (int j = 0; j < tableSize; j++) {
+                    table[j] += (float) ampl * sineTable[sineIndex];
+                    sineIndex += partial;
+                    /* Wrap index at end of table.. */
+                    if (sineIndex >= cycleSize) {
+                        sineIndex -= cycleSize;
+                    }
+                }
+            }
+        }
+
+        /* Normalize after */
+        for (int i = 1; i < numTables; i++) {
+            normalizeArray(tables[i]);
+        }
+    }
+
+    /**************************************************************************/
+    public static float normalizeArray(float[] fdata) {
+        float max, val, gain;
+        int i;
+
+        // determine maximum value.
+        max = 0.0f;
+        for (i = 0; i < fdata.length; i++) {
+            val = Math.abs(fdata[i]);
+            if (val > max)
+                max = val;
+        }
+        if (max < 0.0000001f)
+            max = 0.0000001f;
+        // scale array
+        gain = 1.0f / max;
+        for (i = 0; i < fdata.length; i++)
+            fdata[i] *= gain;
+        return gain;
+    }
+
+    /*****************************************************************************
+     * When the phaseInc maps to the highest level table, then we start interpolating between the
+     * highest table and the raw sawtooth value (phase). When phaseInc points to highest table:
+     * flevel = NUM_TABLES - 1 = -1 - log2(pInc); log2(pInc) = - NUM_TABLES pInc = 2**(-NUM_TABLES)
+     */
+    private final static double LOWEST_PHASE_INC_INV = (1 << NUM_TABLES);
+
+    /**************************************************************************/
+    /* Phase ranges from -1.0 to +1.0 */
+    public double calculateSawtooth(double currentPhase, double positivePhaseIncrement,
+            double flevel) {
+        float[] tableBase;
+        double val;
+        double hiSam; /* Use when verticalFraction is 1.0 */
+        double loSam; /* Use when verticalFraction is 0.0 */
+        double sam1, sam2;
+
+        /* Use Phase to determine sampleIndex into table. */
+        double findex = ((phaseScalar * currentPhase) + phaseScalar);
+        // findex is > 0 so we do not need to call floor().
+        int sampleIndex = (int) findex;
+        double horizontalFraction = findex - sampleIndex;
+        int tableIndex = (int) flevel;
+
+        if (tableIndex > (NUM_TABLES - 2)) {
+            /*
+             * Just use top table and mix with arithmetic sawtooth if below lowest frequency.
+             * Generate new fraction for interpolating between 0.0 and lowest table frequency.
+             */
+            double fraction = positivePhaseIncrement * LOWEST_PHASE_INC_INV;
+            tableBase = tables[(NUM_TABLES - 1)];
+
+            /* Get adjacent samples. Assume guard point present. */
+            sam1 = tableBase[sampleIndex];
+            sam2 = tableBase[sampleIndex + 1];
+            /* Interpolate between adjacent samples. */
+            loSam = sam1 + (horizontalFraction * (sam2 - sam1));
+
+            /* Use arithmetic version for low frequencies. */
+            /* fraction is 0.0 at 0 Hz */
+            val = currentPhase + (fraction * (loSam - currentPhase));
+        } else {
+
+            double verticalFraction = flevel - tableIndex;
+
+            if (tableIndex < 0) {
+                if (tableIndex < -1) // above Nyquist!
+                {
+                    val = 0.0;
+                } else {
+                    /*
+                     * At top of supported range, interpolate between 0.0 and first partial.
+                     */
+                    tableBase = tables[0]; /* Sine wave table. */
+
+                    /* Get adjacent samples. Assume guard point present. */
+                    sam1 = tableBase[sampleIndex];
+                    sam2 = tableBase[sampleIndex + 1];
+
+                    /* Interpolate between adjacent samples. */
+                    hiSam = sam1 + (horizontalFraction * (sam2 - sam1));
+                    /* loSam = 0.0 */
+                    // verticalFraction is 0.0 at Nyquist
+                    val = verticalFraction * hiSam;
+                }
+            } else {
+                /*
+                 * Interpolate between adjacent levels to prevent harmonics from popping.
+                 */
+                tableBase = tables[tableIndex + 1];
+
+                /* Get adjacent samples. Assume guard point present. */
+                sam1 = tableBase[sampleIndex];
+                sam2 = tableBase[sampleIndex + 1];
+
+                /* Interpolate between adjacent samples. */
+                hiSam = sam1 + (horizontalFraction * (sam2 - sam1));
+
+                /* Get adjacent samples. Assume guard point present. */
+                tableBase = tables[tableIndex];
+                sam1 = tableBase[sampleIndex];
+                sam2 = tableBase[sampleIndex + 1];
+
+                /* Interpolate between adjacent samples. */
+                loSam = sam1 + (horizontalFraction * (sam2 - sam1));
+
+                val = loSam + (verticalFraction * (hiSam - loSam));
+            }
+        }
+        return val;
+    }
+
+    public double convertPhaseIncrementToLevel(double positivePhaseIncrement) {
+        if (positivePhaseIncrement < 1.0e-30) {
+            positivePhaseIncrement = 1.0e-30;
+        }
+        return -1.0 - (Math.log(positivePhaseIncrement) / Math.log(2.0));
+    }
+
+    public static MultiTable getInstance() {
+        return instance;
+    }
+
+}