[wozmaker] / src / dsp.cpp

//
// dsp.cpp: Digital Signal Processing module
//
// Part of the WOZ Maker project
// by James Hammons
// (C) 2018 Underground Software
//

#include "dsp.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "fileio.h"
#include "global.h"


//#define NOISY

// Really, there shouldn't be more than 5 waveforms per track
//static uint32_t total = 0;//trkStrms
//static uint32_t slip[9];
//double ratio[9];
uint32_t sNum[9];
double initSyncTime[9];

static const double SIGMA_THRESHOLD = 2.5;


double StdDeviation(double * a, uint32_t num, double * pMean/*= NULL*/)
{
        double sigma = 0, mean = 0;

        for(uint32_t i=0; i<num; i++)
                mean += a[i];

        mean = mean / (double)num;

        if (pMean)
                *pMean = mean;

        for(uint32_t i=0; i<num; i++)
                sigma += (a[i] - mean) * (a[i] - mean);

        sigma *= 1.0 / (double)num;
        return sqrt(sigma);
}


double StdDeviationWave(uint32_t * sync, uint32_t num, double * m = NULL);
double StdDeviationWave(uint32_t * sync, uint32_t num, double * m/*= NULL*/)
{
        double sigma = 0, mean = 0;
        double a[7];

        for(uint32_t i=0; i<num; i++)
        {
                a[i] = (double)Global::wave[sNum[i]][sync[i]];
                mean += a[i];
        }

        mean = mean / (double)num;

        // Capture the mean, if a pointer was passed in
        if (m)
                *m = mean;

        for(uint32_t i=0; i<num; i++)
                sigma += (a[i] - mean) * (a[i] - mean);

        return sqrt(sigma / (double)num);
}


double Largest(double * a, uint32_t num)
{
        double largest = a[0];

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] > largest)
                        largest = a[i];
        }

        return largest;
}


uint32_t Largest(uint32_t * a, uint32_t num)
{
        uint32_t largest = a[0];

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] > largest)
                        largest = a[i];
        }

        return largest;
}


uint32_t LargestIndex(double * a, uint32_t num)
{
        double largest = a[0];
        uint32_t index = 0;

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] > largest)
                {
                        largest = a[i];
                        index = i;
                }
        }

        return index;
}


uint32_t LargestIndex(uint32_t * a, uint32_t num)
{
        uint32_t largest = a[0];
        uint32_t index = 0;

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] > largest)
                {
                        largest = a[i];
                        index = i;
                }
        }

        return index;
}


double Smallest(double * a, uint32_t num)
{
        double smallest = a[0];

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] < smallest)
                        smallest = a[i];
        }

        return smallest;
}


uint32_t Smallest(uint32_t * a, uint32_t num)
{
        uint32_t smallest = a[0];

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] < smallest)
                        smallest = a[i];
        }

        return smallest;
}


uint32_t SmallestIndex(double * a, uint32_t num)
{
        double smallest = a[0];
        uint32_t index = 0;

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] < smallest)
                {
                        smallest = a[i];
                        index = i;
                }
        }

        return index;
}


uint32_t SmallestIndex(uint32_t * a, uint32_t num)
{
        uint32_t smallest = a[0];
        uint32_t index = 0;

        for(uint32_t i=1; i<num; i++)
        {
                if (a[i] < smallest)
                {
                        smallest = a[i];
                        index = i;
                }
        }

        return index;
}


//
// This is where we attempt to use the stream data to collapse them down to a single track.  Taken as wholes, each individual track will drift out of sync with respect to the others, and fairly quickly as the drive mechanism does not read consistently at a set RPM as you might expect, but instead it fluctuates all over the place--which throws a big ol' monkey wrench into the works.  (Yes, the overall *average* stays around a set RPM, but locally--nope.)
//
// The big problem is that the disk spins via a motor which isn't perfect in its rotation that is connected to the hub that spins the disk via a belt, which further adds more uncertainty to the speed of the disk as it rotates.  The end result is that there is no consistency vis-a-vis the speed of the disk, and so the results of one read to the next can and will vary wildly.  It's a miracle that these things are readable at all.  O_o
//
// So how do we cope with it?  I'm sure there's a well known paper or set of papers which addresses noise rejection from multiple signals read from a single source, but I have yet to find them.  So what we have below is a set of ad-hoc assumptions and heuristics that attempt to deal with this problem.
//
// The basic assumption is that from any given pulse to the next one, the drift will tend to be small and thus as we move from pulse to pulse, we can stay more or less in sync by resetting the base from which we calculate the next time interval to the previous step.  This works remarkably well until we run into a long period of zeroes, which seems to throw off the Applesauce hardware which read these streams in the first place; or when some spurious one pulses pop up in one stream that doesn't in one or more of the others.
//
#if 0
uint32_t Synthesize(uint32_t * sync, uint32_t num, uint32_t * synthWave, float * swAmplitude)
{
        uint32_t swLen = 0;
        uint32_t syncPt[9], syncSave[9], nrSync[9];
        double a[9], currentRow[9], nextRow[9];
        double time[9][32], drift[9];
        uint32_t timePtr = 0;

        memcpy(syncPt, sync, sizeof(syncPt));

        while (true)
        {
                for(uint32_t i=0; i<num; i++)
                {
                        if (syncPt[i] >= Global::streamLen[i])
                                return swLen;
                }

                // Save these *before* we step forward
                memcpy(syncSave, syncPt, sizeof(syncSave));
                StepForward(syncPt, a, num);
                double mean;
                double sigma = StdDeviation(a, num, &mean);

                if (sigma <= SIGMA_THRESHOLD)
                {
                        for(uint32_t i=0; i<num; i++)
                                time[i][timePtr] = a[i];

                        timePtr = (timePtr + 1) % 32;

                        synthWave[swLen] = (uint32_t)(mean + 0.5);
                        swAmplitude[swLen] = 1.0f;
                        swLen++;
                }
                else
                {
                        // Calculate drift over the last 32 time steps...
                        for(uint32_t i=0; i<num; i++)
                        {
                                drift[i] = 0;

                                for(uint32_t j=0; j<32; j++)
                                        drift[i] += time[i][j];
                        }

                        double driftSigma = StdDeviation(drift, num);
                        memcpy(nrSync, syncPt, sizeof(nrSync));
                        StepForward(nrSync, nextRow, num);
                        double nextSigma = StdDeviation(nextRow, num);
                        uint32_t lookahead = Lookahead(nrSync, num, +1);

                        // The basic idea is to find the extent of this anomalous region and then use what's there as an amplitude by summing the signals where they correspond.
//                      memcpy(syncSave, syncPt, sizeof(syncSave));
                        bool pulled, doSync = true;
                        uint32_t syncSteps = 0;

                        // Basically we move forward thru the area until either 1) the sigma is OK, or 2) we've pulled everything as far as it can go and thus this becomes our new local start. 3) will go in at some point which is, do a lookahead and see what the sigma looks like; if it's good and the current one is stinko, (not *too* stinko) then we've found a new local start point, otherwise keep going.
/*
Note that 3) works OK to get out of sticky situations, but sometimes the signals are fairly noisy and you will end up in situations where the lookahead is zero and the next sigma isn't great (& the current isn't any good either), but using that as the new local start is OK and will permit walking forward/backward through the stream.  Some of the time you will find the average of the current sigma and the following will be below the threshold with a good lookahead value (10, 100, 1000), this is also a good time to set the new local start.
*/
#ifdef NOISY
printf("a: ");
for(uint32_t i=0; i<num; i++)
        printf("[%.2lf]", a[i]);
printf(" => %.2lf {drift=%.2lf, la=%d}\n", sigma, driftSigma, lookahead);
#endif
                        // New heuristic: if the sigma of the current row is less than the sigma of the drift over the last 32 steps, we consider that to be "in sync"; otherwise, do the usual stuff
/*
This works until you hit this:
a: [1167.00][1174.00][1169.00] => 2.94
          ~~~~
<1167.00><1174.00><1169.00> ==> 2.94 {drift=0.82}(!!!) (lookahead=2.94, 630)
{12168}-{12166}-{12166}
          ~~~~
<1198.00><1205.00><1200.00> ==> 2.94 {drift=0.82} (lookahead=0.00, 629)
{12169}-{12167}-{12167}
          ~~~~
<1231.00><1238.00><1233.00> ==> 2.94 {drift=0.82} (lookahead=0.00, 628)
----------------
<337.00><334.00><335.00> ==> 1.25 (OK) (lookahead=18.40, 1705)
{1260}-{1259}-{1258}
a: [1111.00][1100.00][1097.00] => 6.02 {drift=7.72, la=0}
a: [42.00][51.00][45.00] => 3.74 {drift=7.72, la=150}
  20.00   ~~~~     0.00
< 62.00>< 51.00>< 45.00> ==> 7.04 (lookahead=13.42, 0)
{2972}-{2970}-{2969}
  ~~~~    20.00   13.00 <-- the more I look at it, the more I think we need to coalesce these into their next steps... (can be done in StepForward/Back())
-------------------
it does OK but fails because we don't check *inside* the loop
<363.00><364.00><364.00> ==> 0.47 (OK) (lookahead=52.33, 1704)
{25947}-{25931}-{25922}
a: [33.00][15.00][47.00] => 13.10 {drift=3.86, la=150}
 1149.00  1158.00   ~~~~
<1182.00><1173.00>< 47.00> ==> 532.94 (lookahead=40.94, 0)
{27657}-{27641}-{27627}
  ~~~~    31.00  1155.00
<1182.00><1204.00><1202.00> ==> 9.93 (lookahead=46.20, 6)
{27657}-{27642}-{27633}
  31.00   ~~~~     0.00
<1213.00><1204.00><1202.00> ==> 4.78 (lookahead=0.00, 149)
{27658}-{27642}-{27633}

*/
                        // We check for large zero areas for the new heuristic
//buuuuuuuut... it fails for the stuff that tags along at the end, so...
/*
After looking at this long and hard, and thinking about the fairly random nature of the disk speed as it turns, I've come to the conclusion that trying to use the local speed drift (whatever that means) is futile.  It's random; it just doesn't correlate, period.  So stuff like the following is pretty much crap.
*/
                        if ((mean > 9.0)//00.0)
                                && ((sigma < StdDeviation(drift, num))
                                        || (((sigma + nextSigma) / 2.0) < SIGMA_THRESHOLD)))
                        {
                                doSync = false;
                        }

                        if (doSync)
                        {
                                do
                                {
                                        syncSteps++;

                                        // If we've been in this loop too long, something's very wrong
                                        if (syncSteps > 20)
                                                return swLen;

                                        // Pull the lesser streams forward to catch up to the largest
                                        uint32_t largestIdx = LargestIdx(a, num);
                                        pulled = false;

                                        for(uint32_t i=0; i<num; i++)
                                        {
                                                if (i == largestIdx)
#ifdef NOISY
                                                {
                                                        printf("  ~~~~  ");
#endif
                                                        continue;
#ifdef NOISY
                                                }
#endif

                                                double ttn = 0;
                                                uint32_t walkahead = syncPt[i];

                                                while ((a[largestIdx] - (a[i] + ttn)) > 7.0)
                                                {
                                                        pulled = true;

                                                        if (walkahead >= Global::streamLen[i])
                                                        {
                                                                // We ran out of road, so we're done for now
                                                                return swLen;
                                                        }
/*
Maybe rename this to StepForwardSngl and the other to StepForwardAll?
because otherwise, this can become murky
[well, should be able to tell from the context.
well, the thing that's murky is the all clears the array while the single adds to the existing.  now *that's* confusing...  :-P]
we can call it StepAndAddForward/Back; that would clear it up.
*/
                                                        StepAndAddForward(walkahead, ttn, i);
                                                }

                                                syncPt[i] = walkahead;
                                                a[i] += ttn;
#ifdef NOISY
                                                if (pulled)
                                                        printf(" %6.2lf ", ttn);
                                                else
                                                        printf("        ");
#endif
                                        }

#ifdef NOISY
                                        printf("\n");

                                        for(uint32_t i=0; i<num; i++)
                                                printf("<%6.2lf>", a[i]);
#endif
                                        // Calculate the standard deviation at this new point
                                        sigma = StdDeviation(a, num, &mean);
#ifdef NOISY
                                        printf(" ==> %.2lf%s", sigma, (sigma <= SIGMA_THRESHOLD ? " (OK)" : (!pulled ? "(STEP)" : "")));

double lookahead[9];
//printf(" {drift=%.2lf}", driftSigma);
uint32_t minSyncIdx = SmallestIdx(syncPt, num);
if (syncPt[minSyncIdx] > 0)
{
        for(uint32_t i=0; i<num; i++)
                lookahead[i] = (double)Global::streamData[i][syncPt[i] + 0] * ratio[i];
        printf(" (lookahead=%.2lf, %d)\n", StdDeviation(lookahead, num), Lookahead(syncPt, num, +1));
}
else
        printf("\n");

for(uint32_t j=0; j<num; j++) { printf("{%d}%s", syncPt[j], (j < num - 1 ? "-" : "")); }
printf("\n");
#endif

                                        // Heuristic: see if we're stuck with a bad sigma and can't pull any more forward; if so, pull forward one step and see if that gets us any traction...
                                        if ((sigma > SIGMA_THRESHOLD) && !pulled)
                                        {
                                                for(uint32_t i=0; i<num; i++)
                                                        StepAndAddForward(syncPt[i], a[i], i);

                                                sigma = StdDeviation(a, num, &mean);
                                                pulled = true;
                                        }
                                }
                                while ((sigma > SIGMA_THRESHOLD) && pulled);
                        }

                        // Now, we've re-established sync, now fill the appropriate spots

/*
So the algorithm to do this proppaly is like so:

Step 1: Find the smallest step forward through the section; set amplitude to 0.2.
Step 2: Loop thru the other streams and see if they correspond.  If so, advance their stream start & add 0.2 to the amplitude for each one seen.
Step 3: Add the mean value of the corresponding values to the synthesized waveform + wave amplitude
Step 4: If all of the stream starts are at the end of the section, we are done.
[maybe.  how do you handle the case where you have streams that look like so:

200
150 50
150 50
175 25
200

.......................................|
.............................|.........|
.............................|.........|
..................................|....|
.......................................|

this will fail to add in the correct amount for the times that had no part in the anomaly.  so to fix this, maybe have a baseline value that bumps upward to show where the "cursor" (so to speak) is, so at the end we can subtract the mean that we got above from this baseline to get the remaining time to synthesize, like so:

synthWave[swLen] = (uint32_t)((mean - baseline) + 0.5)
swAmplitude[swLen] = amplitude;

no need for that.  the times that step directly to the end can safely be ignored.
]

150 -> 150 (ignores 175)
150 goes in, ampl. .4, we now have:

200
50
50
175 25
200

only one is left that we can grab:

175

last pulse was at 150 (40%), so we need to subtract that out

-> 25

so maybe we'd do like so:
-1) create a list containing the start values for the start of the anomalous section and replace them as the counters advance.
0) are streams at the end? if so, stuff in the final value to the synthesized wave & stop, otherwise:
1) find the smallest in the current row and use that as a basis.
2) find streams that correlate, if any
3) advance the streams that correlate & subtract out the minimum from the other rows
4) go to 0)

200
150 50
150 50
175 25
200

smallest is 150: pick 150 -> picks up 150
advance the 150 rows & subtract 150 from the rest:

50
50
50
25 25
50

pick smallest: 25 -> picks up nothing
advance the 25 & subtract it from the rest:

25
25
25
25
25

We've reached the end of the section, add the final amount

--------------------------------------------------------------------------------
syncSave before: [2->4][2->6][2->4][2->6][2->6]
syncSave after: [3->4][3->6][3->4][3->6][3->6]

currentRow: [120.89][49.97][46.98][49.98][48.98] => 28.78

2: [3->4][3->6][3->4][3->6][3->6]

a[0]=46.98: [120.89] ==> 36.95
        [49.97]{17.00} ==> 1.50
        <106.00>[49.97][49.98]{28.00} ==> 1.41
        [49.97][49.98][48.98]{26.00} ==> 1.22

3: [3->4][4->6][4->4][4->6][4->6]

a[0]=17.00: [73.91] ==> 28.45
        <56.00>[28.00]{45.00} ==> 5.50
        [28.00][26.00]{47.00} ==> 4.78

4: [3->4][5->6][4->4][5->6][5->6]

a[0]=45.00: [56.91] ==> 5.95
        [56.00]{32.00} ==> 5.50
        <32.00>[56.00][47.00]{32.00} ==> 4.78

5: [3->4][6->6][4->4][6->6][6->6]

a[0]=11.91: <32.00>

6: [4->4][6->6][4->4][6->6][6->6]

 2->4    2->6    2->4    2->6    2->6
[120.89][ 49.97][ 46.98][ 49.98][ 48.98] => 28.78
( 73.91)  17.00  106.00   28.00   26.00
( 56.91)  56.00           45.00   47.00
  32.00   32.00           32.00   32.00

*/
/*printf("syncSave before: ");
for(uint32_t i=0; i<num; i++)
        printf("[%d->%d]", syncSave[i], syncPt[i]);
printf("\n");//*/

                        // Collect the current row
                        StepForward(syncSave, currentRow, num);
                        bool done2;
                        float ampStep = 1.0f / (float)num;

/*printf("syncSave after: ");
for(uint32_t i=0; i<num; i++)
        printf("[%d->%d]", syncSave[i], syncPt[i]);
printf("\n");
printf("currentRow: ");
for(uint32_t i=0; i<num; i++)
        printf("[%.2lf]", currentRow[i]);
printf(" => %.2lf\n", StdDeviation(currentRow, num));//*/

                        done2 = true;

//printf("%d: ", swLen);
                        for(uint32_t i=0; i<num; i++)
                        {
//printf("[%d->%d]", syncSave[i], syncPt[i]);
                                if (syncSave[i] != syncPt[i])
                                {
                                        done2 = false;
//                                      break;
                                }
                        }
//printf("\n");

                        while (!done2)
                        {
                                float amplitude = ampStep;
                                uint32_t smallestIdx = SmallestIdx(currentRow, num);
                                uint32_t cnt = 1;
                                a[0] = currentRow[smallestIdx];
//printf("a[0]=%.2lf: ", a[0]);

                                for(uint32_t i=0; i<num; i++)
                                {
                                        // Make sure we haven't gone beyond the end of any of the streams we're looking at
                                        if (syncSave[i] >= Global::streamLen[i])
                                                return swLen;

                                        if (i == smallestIdx)
                                        {
                                                currentRow[i] = 0;
                                                StepAndAddForward(syncSave[i], currentRow[i], i);
//printf("<%.2lf>", currentRow[i]);

                                                // Make sure no spurious signals are looked at
                                                if ((syncSave[i] < syncPt[i]) && (currentRow[i] < 24.0))
                                                        StepAndAddForward(syncSave[i], currentRow[i], i);

                                                continue;
                                        }

                                        double crSave = currentRow[i];
                                        currentRow[i] -= a[0];

                                        // Don't bother with stuff that's out of range!  :-P
                                        if (syncSave[i] == syncPt[i])
                                                continue;

                                        // Try to find correlates
                                        a[cnt] = crSave;//currentRow[i];
                                        cnt++;
                                        sigma = StdDeviation(a, cnt, &mean);
//for(uint32_t j=1; j<cnt; j++)
//      printf("[%.2lf]", a[j]);

//relaxed sigma...
                                        if (sigma <= 5.5)//SIGMA_THRESHOLD)
                                        {
                                                currentRow[i] = 0;
                                                StepAndAddForward(syncSave[i], currentRow[i], i);
//printf("{%.2lf}", currentRow[i]);
                                                amplitude += ampStep;

                                                // Make sure no spurious signals are looked at
                                                if ((syncSave[i] < syncPt[i]) && (currentRow[i] < 24.0))
                                                        StepAndAddForward(syncSave[i], currentRow[i], i);
                                        }
                                        else
                                        {
                                                cnt--;
                                        }
//printf(" ==> %.2lf\n\t", sigma);
                                }

                                sigma = StdDeviation(a, cnt, &mean);
                                synthWave[swLen] = (uint32_t)(mean + 0.5);
                                swAmplitude[swLen] = amplitude;
                                swLen++;
                                done2 = true;

//printf("%d: ", swLen);
                                for(uint32_t i=0; i<num; i++)
                                {
//printf("[%d->%d]", syncSave[i], syncPt[i]);
                                        if (syncSave[i] != syncPt[i])
                                        {
                                                done2 = false;
//                                              break;
                                        }
                                }
//printf("\n");
                        }

                        sigma = StdDeviation(currentRow, num, &mean);
                        synthWave[swLen] = mean;
                        swAmplitude[swLen] = 1.0f;
                        swLen++;
                }
        }

        return swLen;
}


//
// Attempt to find the sync point for the passed in stream #s.  If found, set
// sync1 & sync2 to the spots where sync was found (typically, one or the other
// will be zero) and return true, otherwise return false.
//
bool FindSync(uint32_t stream1, uint32_t stream2, uint32_t & sync1, uint32_t & sync2)
{
        sync1 = sync2 = 0;

        // We *should* be able to synchronize within around 80 fluxes or so.  If not, then there's bigger problems.
        for(uint32_t offset=0; offset<80; offset++)
        {
                uint32_t strmLen = Uint32LE(Global::stream[stream1]->dataLength);
                bool found = true;

                // How many bits should we consider as synchronized? All of them?
                for(uint32_t i=1, i2=1; i<strmLen; i++, i2++)
                {
                        double time1 = (double)Global::stream[stream1]->data[i];
                        double time2 = (double)Global::stream[stream2]->data[i2 + offset];

                        while (Global::stream[stream1]->data[i] == 0xFF)
                        {
                                i++;
                                time1 += (double)Global::stream[stream1]->data[i];
                        }

                        while (Global::stream[stream2]->data[i2] == 0xFF)
                        {
                                i2++;
                                time2 += (double)Global::stream[stream2]->data[i2 + offset];
                        }

                        double mean = (time1 + time2) * 0.5;
                        double sigma = sqrt(0.5 * (((time1 - mean) * (time1 - mean))
                                + ((time2 - mean) * (time2 - mean))));

                        if (sigma > 2.50)
                        {
                                // How many bits is enough to say we found it?
/*
N.B.: This (1,000) is fine for Apple Panic (which has shitty captures with junk in between sectors), but not for Championship Lode Runner which doesn't have any sync bytes for a loooooooooooong time
*/
                                if (i < 10000)
                                        found = false;

                                break;
                        }
                }

                if (found)
                {
                        sync2 = offset;
                        return true;
                }
        }

        // Since we didn't find sync moving stream 2 forward, let's try stream 1
        for(uint32_t offset=0; offset<80; offset++)
        {
                uint32_t strmLen = Uint32LE(Global::stream[stream2]->dataLength);
                bool found = true;

                // How many bits should we consider as synchronized? All of them?
                for(uint32_t i=1, i2=1; i<strmLen; i++, i2++)
                {
                        double time1 = (double)Global::stream[stream2]->data[i];
                        double time2 = (double)Global::stream[stream1]->data[i2 + offset];

                        while (Global::stream[stream2]->data[i] == 0xFF)
                        {
                                i++;
                                time1 += (double)Global::stream[stream2]->data[i];
                        }

                        while (Global::stream[stream1]->data[i2] == 0xFF)
                        {
                                i2++;
                                time2 += (double)Global::stream[stream1]->data[i2 + offset];
                        }

                        double mean = (time1 + time2) * 0.5;
                        double sigma = sqrt(0.5 * (((time1 - mean) * (time1 - mean))
                                + ((time2 - mean) * (time2 - mean))));

                        if (sigma > 2.50)
                        {
                                // How many bits is enough to say we found it?
                                if (i < 10000)
                                        found = false;

                                break;
                        }
                }

                if (found)
                {
                        sync1 = offset;
                        return true;
                }
        }

        return false;
}


void FindSyncForStreams(uint32_t * syncPt, uint32_t num)
{
//      uint32_t syncPt[10] = { 0 };
        uint32_t base = 0;
        syncPt[0] = 0;

        for(uint32_t i=1; i<num; i++)
        {
                if (FindSync(sNum[base], sNum[i], syncPt[base], syncPt[i]) == true)
                {
                        // Check to see which stream has to move.  If it's not the current base, then we don't have to do anything but check the next one; otherwise, we have to shift everything done prior by the amount of the shift and set the new base to be the one that didn't move.
                        if (syncPt[base] != 0)
                        {
                                for(uint32_t j=0; j<i; j++)
                                {
                                        if (j != base)
                                                syncPt[j] += syncPt[base];
                                }

                                base = i;
                        }
                }
                else
                {
                        // should we do this? or just set a flag that we couldn't sync this one?
                        syncPt[i] = -1;
//                      return;
                }
        }

        // If there was no sync for any of the other tracks, signal this by counting the number of non-sync tracks and checking against the total # of tracks - 1.
        uint32_t nonSync = 0;

        for(uint32_t i=1; i<num; i++)
        {
                if (syncPt[i] == (uint32_t)-1)
                        nonSync++;
        }

        // Signal that sync failed
        if (nonSync == (num - 1))
                syncPt[0] = -1;
/*
v v
0 0 0 0 0  base = 0, i = 1

0  1
0, 10

v    v
0 10 0 0 0  base = 0, i = 2
5 10 0 0 0

0  2
5, 0

     v v
5 15 0 0 0  base = 2, i = 3
5 15 0 7 0

2  3
0, 7

     v   v
5 15 0 7 0  base = 2, i = 4
5 15 8 7 0

2  4
8, 0

13 23 8 15 0  base = 4

        if (FindSync(sNum[0], sNum[1], syncPt[0], syncPt[1]) == false)
        {
                // give up, we couldn't sync stream 0 & 1...  :-/
        }

        // check to see which one had to move, and use the other to check against the rest...
        uint32_t base = (syncPt[0] == 0 ? 0 : 1);

        if (FindSync(sNum[base], sNum[2], syncPt[base], syncPt[2]) == false)
        {
                // give up, we couldn't sync stream 0 & 1...  :-/
        }

        // check again, to see which one moved and the other...
        uint32_t base2 = (syncPt[base] == 0 ? base, 2);

        // & so on...


*/
}


void FindInitialSyncForStreams(uint32_t * sync, uint32_t num)
{
/*
What we need to do is find the FF. or FF.. or FE..s, then find the first non-zero trailing byte after that to be our sync point.  If we don't, we run the risk of getting bad starts because some bitstreams might have some sync starts that start with the wrong number of 1 bits, and thus walking back will fail.
*/
        // Let's try aligning samples to each other by counting to a short stream
        // of FF.(.)s...
        uint64_t bitPat0 = 0;

        for(uint32_t i=0; i<num; i++)
        {
                sync[i] = 0;
                uint64_t bitPat = 0;
                initSyncTime[i] = 0;

                // We'll use a naïve window to scan for the bit patterns
                for(int32_t j=0; j<Uint32LE(Global::stream[sNum[i]]->dataLength); j++)
                {
                        uint32_t timeToNext = Global::streamData[i][j];

                        while (Global::streamData[i][j] == 0xFF)
                        {
                                j++;
                                timeToNext += Global::streamData[i][j];
                        }

                        initSyncTime[i] += (double)timeToNext;

                        if ((timeToNext >= 24) && (timeToNext <= 49))
                                bitPat = (bitPat << 1) | 1;
                        else if ((timeToNext >= 50) && (timeToNext <= 80))
                                bitPat = (bitPat << 2) | 1;
                        else if ((timeToNext >= 81) && (timeToNext <= 112))
                                bitPat = (bitPat << 3) | 1;
                        else
                        {
                                while (timeToNext > 32)
                                {
                                        bitPat <<= 1;
                                        timeToNext -= 32;
                                }
                        }

                        if (i == 0)
                                bitPat0 = bitPat;

                        // Search for first instance of 3 FF. or 3 FF.. in bitstream
                        // Might need to look for FE.. as well
                        //    0 1111 1111 0111 1111 1011 1111 1101
                        // 0111 1111 1001 1111 1110 0111 1111 1001
//                      if (((bitPat & 0x1FFFFFFF) == 0x0FF7FBFD)
//                              || (bitPat & 0xFFFFFFFF) == 0x7F9FE7F9)
                        // Need to see at least two regular bytes after as well:
                        //    0 1111 1111 0111 1111 1011 1111 1101 xxxx xxx1 xxxx xxx1
                        // 0111 1111 1001 1111 1110 0111 1111 1001 xxxx xxx1 xxxx xxx1
                        if ((bitPat & 0x1FFFFFFF0101) == 0x0FF7FBFD0101)
                        {
                                if ((i != 0) && (bitPat != bitPat0))
                                        continue;

                                sync[i] = j - (j - 25 >= 0 ? 25 : 0);
                                break;
                        }
                        else if ((bitPat & 0xFFFFFFFF0101) == 0x7F9FE7F90101)
                        {
                                if ((i != 0) && (bitPat != bitPat0))
                                        continue;

                                sync[i] = j - (j - 25 >= 0 ? 25 : 0);
                                break;
                        }
                }
        }
}


//
// Walk through the streams backwards until we hit the zero point on one of them.  This assumes the sync point we found initially is a good one (it might not be!).
//
void BacktrackInitialSync(uint32_t * syncPt, uint32_t num)
{
//      uint32_t count = 0;
        double a[10];
//      uint32_t syncSave[10];

        // Save the original sync points (do we need to?)
//      for(uint32_t i=0; i<num; i++)
//              syncSave[i] = syncPt[i];

        uint32_t minSlip = (uint32_t)-1;

        for(uint32_t i=0; i<num; i++)
        {
                if (syncPt[i] < minSlip)
                        minSlip = syncPt[i];
        }

        while (minSlip > 0)
        {
                // Get the time for this step into a[]
                for(uint32_t i=0; i<num; i++)
                {
                        // We scale the time by a factor of the measured RPM to 300 RPM, which should get us closer to a "normal" time for all streams
                        a[i] = (double)Global::stream[sNum[i]]->data[syncPt[i]];// * ratio[i];

                        while ((syncPt[i] > 0)
                                && (Global::stream[sNum[i]]->data[syncPt[i] - 1] == 0xFF))
                        {
                                syncPt[i]--;
                                a[i] += (double)Global::stream[sNum[i]]->data[syncPt[i]];// * ratio[i];
                        }
                }

                double sigma = StdDeviation(a, num);
//              count++;

                if (sigma > SIGMA_THRESHOLD)
                {
                        printf("\n");

                        for(uint32_t i=0; i<num; i++)
//                              printf("<%6.2lf>", (double)Global::stream[sNum[i]]->data[syncPt[i]]);// * ratio[j]);
                                printf("<%6.2lf>", a[i]);// * ratio[j]);

                        printf(" --> %.2lf\n", sigma);

                        uint32_t syncCount = 0;

                        // Save sync origins for later analysis
//                      for(uint32_t j=0; j<num; j++)
//                              syncSave[j] = syncPt[j];

                        // It may turn out that this is too large too...
                        while (sigma > 3.30)//SIGMA_THRESHOLD)
                        {
                                // Deal with the problem.  The heuristic is find the largest value in the set, then advance the other streams by one step until they are "close" to the largest.  If the std. deviation is too large, try again; eventually it should be within good limits.
                                uint32_t largestIdx = LargestIdx(a, num);

                                for(uint32_t i=0; i<num; i++)
                                {
                                        if (i == largestIdx)
                                        {
                                                printf("  ~~~~  ");
                                                continue;
                                        }

                                        int32_t walkBack = syncPt[i] - 1;
                                        double ttn = (double)Global::stream[sNum[i]]->data[walkBack];

                                        while ((walkBack > 0)
                                                && (Global::stream[sNum[i]]->data[walkBack - 1] == 0xFF))
                                        {
                                                walkBack--;
                                                ttn += (double)Global::stream[sNum[i]]->data[walkBack];// * ratio[i];
                                        }

                                        // If the difference between the highest and this is over 5, add the next time interval and let the sync counter, uh, slip (backwards!) :-)
                                        if ((a[largestIdx] - a[i]) > 5.0)
                                        {
                                                syncPt[i] = walkBack;
                                                a[i] += ttn;
                                                printf(" %6.2lf ", ttn);
                                        }
                                        else
                                                printf("        ");
                                }

                                printf("\n");

                                for(uint32_t i=0; i<num; i++)
                                        printf("<%6.2lf>", a[i]);

                                sigma = StdDeviation(a, num);
                                printf(" ==> %.2lf\n", sigma);

                                syncCount++;

                                if (syncCount > 20)
                                {
                                        // Back up the train
                                        for(uint32_t i=0; i<num; i++)
                                        {
                                                syncPt[i] -= (syncPt[i] >= 40 ? 40 : 0);

                                                if ((int32_t)syncPt[i] < 0)
                                                        syncPt[i] = 0;
                                        }

                                        minSlip = 0;
                                        printf("Could not synchronize streams...\n");
                                        break;
                                }
                        }
                }

                // Point as the next time step and see where the minimum sync point is now...
                for(uint32_t i=0; i<num; i++)
                {
                        syncPt[i]--;

                        if (syncPt[i] < minSlip)
                                minSlip = syncPt[i];
                }
        }
}


//
// Walk through the streams backwards until we hit the zero point on one of them.  This assumes the sync point we found initially is a good one (it might not be!). [We've found that sometimes, the initial sync point is *really* bad.]
//
void BacktrackInitialSync2(uint32_t * syncPt, uint32_t num)
{
        double a[10];
        uint32_t syncSave[10];
        uint32_t minSlip = Smallest(syncPt, num);

        while (minSlip > 0)
        {
                // Get the time for this step into a[]
                for(uint32_t i=0; i<num; i++)
                {
                        // We scale the time by a factor of the measured RPM to 300 RPM, which should get us closer to a "normal" time for all streams
                        a[i] = (double)Global::stream[sNum[i]]->data[syncPt[i]];// * ratio[i];

                        while ((syncPt[i] > 0)
                                && (Global::stream[sNum[i]]->data[syncPt[i] - 1] == 0xFF))
                        {
                                syncPt[i]--;
                                a[i] += (double)Global::stream[sNum[i]]->data[syncPt[i]];// * ratio[i];
                        }
                }

                double sigma = StdDeviation(a, num);

                if (sigma > SIGMA_THRESHOLD)
                {
                        uint32_t syncCount = 0;

                        // Save sync origins for later analysis
                        memcpy(syncSave, syncPt, sizeof(syncSave));

                        // It may turn out that this is too large too...
                        while (sigma > SIGMA_THRESHOLD)
                        {
                                // Deal with the problem.  The heuristic is find the largest value in the set, then advance the other streams by one step until they are "close" to the largest.  If the std. deviation is too large, try again; eventually it should be within good limits.
                                uint32_t largestIdx = LargestIdx(a, num);

                                // We now pull all streams "forward" until they are within 7
                                // underneath the largest or past it.
                                for(uint32_t i=0; i<num; i++)
                                {
                                        if (i == largestIdx)
                                                continue;

                                        double ttn = 0;
                                        int32_t walkback = syncPt[i];

                                        while((a[largestIdx] - (a[i] + ttn)) > 7.0)
                                        {
                                                if (walkback <= 0)
                                                {
                                                        memcpy(syncPt, syncSave, sizeof(syncSave));
                                                        printf("Synchronized as far back as possible...\n");
                                                        return;
                                                }

                                                walkback--;
                                                ttn += (double)Global::stream[sNum[i]]->data[walkback];// * ratio[i];

                                                while ((walkback > 0)
                                                        && (Global::stream[sNum[i]]->data[walkback - 1] == 0xFF))
                                                {
                                                        walkback--;
                                                        ttn += (double)Global::stream[sNum[i]]->data[walkback];// * ratio[i];
                                                }
                                        }

                                        syncPt[i] = walkback;
                                        a[i] += ttn;
                                }

                                sigma = StdDeviation(a, num);

                                syncCount++;

                                if (syncCount > 20)
                                {
                                        // Back up the train
                                        for(uint32_t i=0; i<num; i++)
                                        {
                                                syncPt[i] -= (syncPt[i] >= 40 ? 40 : 0);

                                                if ((int32_t)syncPt[i] < 0)
                                                        syncPt[i] = 0;
                                        }

                                        minSlip = 0;
                                        printf("Could not synchronize streams...\n");
                                        break;
                                }
                        }
                }

                // Point as the next time step and see where the minimum sync point is now...
                for(uint32_t i=0; i<num; i++)
                {
                        syncPt[i]--;

                        if (syncPt[i] < minSlip)
                                minSlip = syncPt[i];
                }
        }
}


//
// Walk through the streams backwards until we hit the zero point on one of them.  This assumes the sync point we found initially is a good one (it might not be!). [We've found that sometimes, the initial sync point is *really* bad.]
//
void BacktrackInitialSync3(uint32_t * sync, uint32_t num)
{
        double a[9];
        uint32_t syncSave[9];

        while (true)
        {
                uint32_t minSlip = Smallest(sync, num);

                if (minSlip == 0)
                        break;

                // Save sync origins for later analysis if necessary
                memcpy(syncSave, sync, sizeof(syncSave));

                // Get the time for this step into a[]
                StepBack(sync, a, num);
                double sigma = StdDeviation(a, num);

                if (sigma > SIGMA_THRESHOLD)
                {
                        uint32_t syncCount = 0;
                        bool pulled;

                        do
                        {
                                syncCount++;

                                if (syncCount > 20)
                                {
                                        // Back up the train
                                        for(uint32_t i=0; i<num; i++)
                                        {
                                                sync[i] -= (sync[i] >= 40 ? 40 : 0);

                                                if ((int32_t)sync[i] < 0)
                                                        sync[i] = 0;
                                        }

                                        printf("Could not synchronize streams...\n");
                                        return;
                                }

                                // Deal with the problem.  The heuristic is find the largest value in the set, then advance the other streams by one step until they are "close" to the largest.  If the std. deviation is too large, try again; eventually it should be within good limits.
                                uint32_t largestIdx = LargestIdx(a, num);
                                pulled = false;

                                // We now pull all streams "forward" until they are within 7
                                // underneath the largest or past it.
                                for(uint32_t i=0; i<num; i++)
                                {
                                        if (i == largestIdx)
                                                continue;

                                        double ttn = 0;
                                        uint32_t walkback = sync[i];

                                        while((a[largestIdx] - (a[i] + ttn)) > 7.0)
                                        {
                                                pulled = true;

                                                if ((int32_t)walkback <= 0)
                                                {
                                                        memcpy(sync, syncSave, sizeof(syncSave));
                                                        printf("Synchronized as far back as possible...\n");
                                                        return;
                                                }

                                                StepAndAddBack(walkback, ttn, i);
                                        }

                                        sync[i] = walkback;
                                        a[i] += ttn;
                                }

                                // Calculate the standard deviation at this new point
                                sigma = StdDeviation(a, num);

                                // Heuristic: see if we're stuck with a bad sigma and can't pull any more forward; if so, pull forward one step and see if that gets us any traction...
                                if ((sigma > SIGMA_THRESHOLD) && !pulled)
                                {
                                        for(uint32_t i=0; i<num; i++)
                                                StepAndAddBack(sync[i], a[i], i);

                                        sigma = StdDeviation(a, num);//, &mean);
                                        pulled = true;
                                }
                        }
                        while ((sigma > SIGMA_THRESHOLD) && pulled);
                }
        }
}


void BacktrackInitialSync4(uint32_t * sync, uint32_t num)
{
        double a[9], a2[9];
        uint32_t syncSave[9];

        while (true)
        {
                uint32_t minSlip = Smallest(sync, num);

                if (minSlip == 0)
                        break;

                // Save sync origins for later analysis if necessary
                memcpy(syncSave, sync, sizeof(syncSave));

                // Get the time for this step into a[]
                StepBack(sync, a, num);
                double sigma = StdDeviation(a, num);

                if (sigma > SIGMA_THRESHOLD)
                {
printf("\n");
for(uint32_t i=0; i<num; i++) { printf("<%6.2lf>", a[i]); }
printf(" --> %.2lf\n", sigma);
                        // Let's see if we can re-sync through this rough patch
                        int result = Resync(sync, a, num, -1);

                        // Well, we failed to find synchronization...  :-(
                        if (result == -1)
                        {
                                memcpy(sync, syncSave, sizeof(syncSave));
                                return;
                        }

                        // OK, we have a good result, but what kind (1=end, 0=sync OK)?
                        if (result == 0)
                                StepBack(sync, a2, num);
                        else if (result == 1)
                        {
for(uint32_t j=0; j<num; j++) { printf("{%d} ", sync[j]); }
printf("\n");
                                memcpy(sync, syncSave, sizeof(syncSave));
for(uint32_t j=0; j<num; j++) { printf("<%d> ", sync[j]); }
printf("\n");
                                StepBack(sync, a2, num);
                                return;
                        }
                }
        }
}


void StepBackUntilBad(uint32_t * sync, double * a, uint32_t num)
{
        while (true)
        {
                // See where the minimum sync point is now...
                uint32_t minSlip = Smallest(sync, num);

                if (minSlip == 0)
                        break;

                // Get the time for this step into a[]
                StepBack(sync, a, num);
                double sigma = StdDeviation(a, num);

                if (sigma > SIGMA_THRESHOLD)
                {
                        printf("\n");

                        for(uint32_t i=0; i<num; i++)
                                printf("<%6.2lf>", a[i]);

                        printf(" --> %.2lf\n", sigma);
                        return;
                }
        }

        printf("---------------------------------------------------\nAT START\n\n");
}


bool StepBackThruBad(uint32_t * sync, double * a, uint32_t num)
{
        // Deal with the problem.  The heuristic is find the largest value in the set, then advance the other streams by one step until they are "close" to the largest.  If the std. deviation is too large, try again; eventually it should be within good limits.
        uint32_t largestIdx = LargestIdx(a, num);

/*
New heuristic to try here: advance each stream until it comes within 7 on the bottom end or passes it altogether; check the sigma then.

2nd heuristic: attempt to pull the next 'num' time all at once and *add* it to the current step, check the sigma to see if it's in tolerance (*and* check the lookahead).
*/
        for(uint32_t i=0; i<num; i++)
        {
                if (i == largestIdx)
                {
                        printf("  ~~~~  ");
                        continue;
                }

                // heuristic 1
                bool seen = false;
                double ttn = 0;
                uint32_t walkback = sync[i];

                // If the difference between the highest and this is over 7, add the next time interval and let the sync counter, uh, slip (backwards!) :-)
                while ((a[largestIdx] - (a[i] + ttn)) > 7.0)
                {
                        seen = true;
                        StepAndAddBack(walkback, ttn, i);
                }

                sync[i] = walkback;
                a[i] += ttn;

                if (seen)
                        printf(" %6.2lf ", ttn);
                else
                        printf("        ");
        }

        printf("\n");

        for(uint32_t i=0; i<num; i++)
                printf("<%6.2lf>", a[i]);

        double sigma = StdDeviation(a, num);
        printf(" ==> %.2lf%s", sigma, (sigma <= SIGMA_THRESHOLD ? " (OK)" : ""));

double lookahead[9];
for(uint32_t i=0; i<num; i++)
        lookahead[i] = (double)Global::streamData[i][sync[i] - 1] * ratio[i];
printf(" (lookahead=%.2lf, %d)\n", StdDeviation(lookahead, num), Lookahead(sync, num, -1));
for(uint32_t j=0; j<num; j++) { printf("{%d}%s", sync[j], (j < num - 1 ? "-" : "")); }
printf("\n");

        if (sigma < SIGMA_THRESHOLD)
        {
                // If we're out of the woods, make sure to step to the next step...
                double b[9];
                StepBack(sync, b, num);
                return false;
        }

        return true;
}


void StepBack(uint32_t * sync, double * a, uint32_t num)
{
        for(uint32_t i=0; i<num; i++)
        {
                int32_t walkback = sync[i] - 1;
                a[i] = (double)Global::streamData[i][walkback] * ratio[i];

                while ((walkback > 0)
                        && (Global::streamData[i][walkback - 1] == 0xFF))
                {
                        walkback--;
                        a[i] += (double)Global::streamData[i][walkback] * ratio[i];
                }

                sync[i] = walkback;
        }
}


void StepAndAddBack(uint32_t * sync, double * a, uint32_t num)
{
        for(uint32_t i=0; i<num; i++)
        {
                int32_t walkback = sync[i] - 1;
                a[i] += (double)Global::streamData[i][walkback] * ratio[i];

                while ((walkback > 0)
                        && (Global::streamData[i][walkback - 1] == 0xFF))
                {
                        walkback--;
                        a[i] += (double)Global::streamData[i][walkback] * ratio[i];
                }

                sync[i] = walkback;
        }
}


void StepAndAddBack(uint32_t & loc, double & a, uint32_t i)
{
        int32_t walkback = loc - 1;
        a += (double)Global::streamData[i][walkback] * ratio[i];

        while ((walkback > 0)
                && (Global::streamData[i][walkback - 1] == 0xFF))
        {
                walkback--;
                a += (double)Global::streamData[i][walkback] * ratio[i];
        }

        loc = (uint32_t)walkback;
}


//
// Move every stream forward by one step
//
void StepForward(uint32_t * sync, double * a, uint32_t num)
{
        for(uint32_t i=0; i<num; i++)
        {
                uint32_t walkahead = sync[i] + 1;
                a[i] = (double)Global::streamData[i][walkahead] * ratio[i];

                while ((walkahead < Global::streamLen[i])
                        && (Global::streamData[i][walkahead] == 0xFF))
                {
                        walkahead++;
                        a[i] += (double)Global::streamData[i][walkahead] * ratio[i];
                }

                sync[i] = walkahead;
        }
}


//
// Move every stream forward by one step
//
void StepAndAddForward(uint32_t * sync, double * a, uint32_t num)
{
        for(uint32_t i=0; i<num; i++)
        {
                uint32_t walkahead = sync[i] + 1;
                a[i] += (double)Global::streamData[i][walkahead] * ratio[i];

                while ((walkahead < Global::streamLen[i])
                        && (Global::streamData[i][walkahead] == 0xFF))
                {
                        walkahead++;
                        a[i] += (double)Global::streamData[i][walkahead] * ratio[i];
                }

                sync[i] = walkahead;
        }
}


//
// Move just *this* stream forward by one step
//
void StepAndAddForward(uint32_t & loc, double & a, uint32_t i)
{
        loc++;
        a += (double)Global::streamData[i][loc] * ratio[i];

        while ((loc < Global::streamLen[i])
                && (Global::streamData[i][loc] == 0xFF))
        {
                loc++;
                a += (double)Global::streamData[i][loc] * ratio[i];
        }
}


//
// Return the number of times with good correlation in the specified direction
// (1 for forward, -1 for backwards)
//
uint32_t Lookahead(uint32_t * syncPt, uint32_t num, int32_t dir /*= 1*/)
{
        uint32_t tSync[10];
        double a[10];
        uint32_t count = 0;

        // Use temporary sync array instead of the one passed in
        memcpy(tSync, syncPt, sizeof(tSync));

        if (dir == -1)
        {
                uint32_t smallestIdx = SmallestIdx(tSync, num);

                while (tSync[smallestIdx] > 0)
                {
                        StepBack(tSync, a, num);
                        double sigma = StdDeviation(a, num);

                        if (sigma > SIGMA_THRESHOLD)
                                break;

                        count++;
                        smallestIdx = SmallestIdx(tSync, num);
                }
        }
        else
        {
                uint32_t largestIdx = LargestIdx(tSync, num);

                while (tSync[largestIdx] < Global::streamLen[largestIdx])
                {
                        StepForward(tSync, a, num);
                        double sigma = StdDeviation(a, num);

                        if (sigma > SIGMA_THRESHOLD)
                                break;

                        count++;
                        largestIdx = LargestIdx(tSync, num);
                }
        }

        return count;
}


//
// Attempt to resynchronize the current streams
// N.B.: When we get here, we're already at the point *after* where it blew up
//       the sigma (the times for that step are passed-in in "a")
//
// Returns -1 on failure, 0 is sync was successful, or 1 if reached end of
// stream before finding sync
//
int Resync(uint32_t * sync, double * a, uint32_t num, int32_t dir/*= 1*/)
{
        uint32_t syncCount = 0;
        bool pulled;
        double sigma = 0;

for(uint32_t j=0; j<num; j++) { printf("{%d} ", sync[j]); }
printf("\n");
        do
        {
                syncCount++;

                if (syncCount > 20)
                {
                        // Back up the train
                        printf("Could not synchronize streams...\n");
                        return -1;
                }

                // Deal with the problem.  The heuristic is find the largest value in the set, then advance the other streams by one step until they are "close" to the largest.  If the std. deviation is too large, try again; eventually it should be within good limits.
                uint32_t largestIdx = LargestIdx(a, num);
                pulled = false;

                // We now pull all streams "forward" until they are within 7
                // underneath the largest or past it.
                for(uint32_t i=0; i<num; i++)
                {
                        if (i == largestIdx)
                        {
                                printf("  ~~~~  ");
                                continue;
                        }

                        double ttn = 0;
                        uint32_t walkback = sync[i];

                        while ((a[largestIdx] - (a[i] + ttn)) > 7.0)
                        {
                                pulled = true;

                                if (((dir == 1) && (walkback >= Global::streamLen[i]))
                                        || ((dir == -1) && ((int32_t)walkback <= 0)))
                                {
//                                      memcpy(sync, syncSave, sizeof(syncSave));
                                        printf("\nSynchronized as far back as possible...\n");
for(uint32_t j=0; j<num; j++) { printf("{%d} ", sync[j]); }
printf("\n");
                                        return 1;
                                }

                                if (dir == 1)
                                        StepAndAddForward(walkback, ttn, i);
                                else
                                        StepAndAddBack(walkback, ttn, i);
                        }

                        sync[i] = walkback;
                        a[i] += ttn;

                        if (pulled)
                                printf(" %6.2lf ", ttn);
                        else
                                printf("        ");
                }

                printf("\n");

                for(uint32_t i=0; i<num; i++)
                        printf("<%6.2lf>", a[i]);

                // Calculate the standard deviation at this new point
                sigma = StdDeviation(a, num);
printf(" ==> %.2lf%s", sigma, (sigma <= SIGMA_THRESHOLD ? " (OK)" : (!pulled ? "(STEP)" : "")));

double lookahead[9];
uint32_t minSyncIdx = SmallestIdx(sync, num);
if (sync[minSyncIdx] > 0 && dir == -1)
{
for(uint32_t i=0; i<num; i++)
{       lookahead[i] = (double)Global::streamData[i][sync[i] + dir] * ratio[i]; }
printf(" (lookahead=%.2lf, %d)\n", StdDeviation(lookahead, num), Lookahead(sync, num, dir));
}
                // Heuristic: see if we're stuck with a bad sigma and can't pull any more forward; if so, pull forward one step and see if that gets us any traction...
                if ((sigma > SIGMA_THRESHOLD) && !pulled)
                {
                        if (dir == 1)
                                StepAndAddForward(sync, a, num);
                        else
                                StepAndAddBack(sync, a, num);

                        sigma = StdDeviation(a, num);//, &mean);
                        pulled = true;
                }
        }
        while ((sigma > SIGMA_THRESHOLD) && pulled);

//      double a2[9];
//      StepBack(sync, a2, num);

        return 0;
}

#endif


//
// Compare two arrays to see if they are equal
//
bool Equal(uint32_t * a1, uint32_t * a2, uint32_t num)
{
        for(uint32_t i=0; i<num; i++)
        {
                if (a1[i] != a2[i])
                        return false;
        }

        return true;
}


uint32_t LookaheadWave(uint32_t * passedInSync, uint32_t num, int8_t dir = 1);
uint32_t LookaheadWave(uint32_t * passedInSync, uint32_t num, int8_t dir/*= 1*/)
{
        uint32_t sync[6];
        uint32_t count = 0;
        memcpy(sync, passedInSync, sizeof(sync));

        while (true)
        {
                if (StdDeviationWave(sync, num) > 2.50)
                        break;

                for(uint32_t i=0; i<num; i++)
                {
                        if (((dir == 1) && (sync[i] >= Global::waveLen[sNum[i]]))
                                || ((dir == -1) && (sync[i] == 0)))
                                return count;

                        sync[i] += dir;
                }

                count++;
        }

        return count;
}


uint32_t FindSyncBetweenFirstTwo(uint32_t * sync)
{
        uint32_t bestSync[2], bestLA = 0, syncSave = sync[0];

        // Move the first stream first...
        for(uint32_t i=0; i<200; i++)
        {
                uint32_t la = LookaheadWave(sync, 2);

                // See if the sync we've found is the best so far...
                if (la > bestLA)
                {
                        bestLA = la;
                        memcpy(bestSync, sync, sizeof(bestSync));
                }

                sync[0]++;
        }

        // Otherwise, reset the first stream and move the 2nd
        sync[0] = syncSave;

        for(uint32_t i=0; i<200; i++)
        {
                // We can bump this first since we already checked the [0][0] case above
                sync[1]++;
                uint32_t la = LookaheadWave(sync, 2);

                // See if the sync we've found is the best so far...
                if (la > bestLA)
                {
                        bestLA = la;
                        memcpy(bestSync, sync, sizeof(bestSync));
                }
        }

        // Set sync to the best we found & return the best lookahead count
        memcpy(sync, bestSync, sizeof(bestSync));

        return bestLA;
}


void FindInitialSyncForStreams2(uint32_t * sync, uint32_t num)
{
        bool lookAt[6] = { true, true, false, false, false, false };

//      while (FindSyncBetweenFirstTwo(sync) == false)
        while (FindSyncBetweenFirstTwo(sync) < 1000)
        {
                for(uint32_t i=0; i<num; i++)
                {
                        sync[i] += 200;

                        if (sync[i] >= Global::waveLen[sNum[i]])
                        {
                                printf("Could not find initial sync!!\n");
                                return;// -1;
                        }
                }
        }

        uint32_t lookahead = LookaheadWave(sync, 2);

        if (lookAt[0])
                printf("Found sync for streams 1 & 2 at %d and %d (la=%d)\n", sync[0], sync[1], lookahead);

        // Add in the others, one at a time and try to synchronize them
        uint32_t bestSync[6];

        for(uint32_t i=2; i<num; i++)
        {
                lookAt[i] = true;
                uint32_t bestLA = 0;

                for(uint32_t j=0; j<200; j++)
                {
                        uint32_t syncSave = sync[i];

                        // And check again for a match
                        for(uint32_t k=0; k<200; k++)
                        {
                                uint32_t la = LookaheadWave(sync, i + 1);

                                if (la > bestLA)
                                {
                                        bestLA = la;
                                        memcpy(bestSync, sync, sizeof(bestSync));
                                }

                                sync[i]++;
                        }

                        sync[i] = syncSave;

                        // Kick synchronized streams forward to check the rest...
                        for(uint32_t k=0; k<i; k++)
                                sync[k]++;
                }

                memcpy(sync, bestSync, sizeof(bestSync));
        }

        for(uint32_t i=0; i<num; i++)
                printf("Sync for stream %d: %d (%s)\n", i, sync[i], (lookAt[i] ? "yes" : "no"));
        printf("Lookahead: %d\n", LookaheadWave(sync, num));
}


#if 1
uint32_t Synthesize2(uint32_t * sync, uint32_t num, uint32_t * synthWave, float * swAmplitude)
{
        uint32_t swLen = 0;
        uint32_t syncPt[7], syncSave[7];
        double a[7], b[7], currentRow[7];
        Global::bloops = 0;

        memcpy(syncPt, sync, sizeof(syncPt));

        while (true)
        {
                for(uint32_t i=0; i<num; i++)
                {
                        if (syncPt[i] >= Global::streamLen[i])
                                return swLen;
                }

                // Save these *before* we step forward
                memcpy(syncSave, syncPt, sizeof(syncSave));
//              StepForward(syncPt, a, num);
                for(uint32_t i=0; i<num; i++)
                {
                        a[i] = (double)Global::wave[sNum[i]][syncPt[i]];
                        syncPt[i]++;
                }

                double mean;
                double sigma = StdDeviation(a, num, &mean);

                if (sigma <= SIGMA_THRESHOLD)
                {
                        synthWave[swLen] = (uint32_t)(mean + 0.5);
                        swAmplitude[swLen] = 1.0f;
                        swLen++;
                }
                else
                {
                        Global::bloops++;
                        // The basic idea is to find the extent of this anomalous region and then use what's there as an amplitude by summing the signals where they correspond.
                        uint32_t syncSteps = 0;
                        uint32_t lookahead = LookaheadWave(syncPt, num);

                        // Basically we move forward thru the area until either 1) the sigma is OK, or 2) we've pulled everything as far as it can go and thus this becomes our new local start. 3) will go in at some point which is, do a lookahead and see what the sigma looks like; if it's good and the current one is stinko, (not *too* stinko) then we've found a new local start point, otherwise keep going.
/*
Note that 3) works OK to get out of sticky situations, but sometimes the signals are fairly noisy and you will end up in situations where the lookahead is zero and the next sigma isn't great (& the current isn't any good either), but using that as the new local start is OK and will permit walking forward/backward through the stream.  Some of the time you will find the average of the current sigma and the following will be below the threshold with a good lookahead value (10, 100, 1000), this is also a good time to set the new local start.
*/
#ifdef NOISY
printf("a: ");
for(uint32_t i=0; i<num; i++)
        printf("[%.2lf]", a[i]);
//printf(" => %.2lf {drift=%.2lf, la=%d}\n", sigma, driftSigma, lookahead);
printf(" => %.2lf {la=%d}\n", sigma, lookahead);

for(uint32_t j=0; j<num; j++) { printf("{%d}%s", syncPt[j], (j < num - 1 ? "-" : "")); }
printf("\n");//*/
#endif
                        bool keepWorking = true;

                        do
                        {
                                syncSteps++;

                                // If we've been in this loop too long, something's very wrong
                                if (syncSteps > 20)
                                        return swLen;

                                // Pull the lesser streams forward to catch up to the largest
                                uint32_t largestIdx = LargestIndex(a, num);
                                bool pulled = false;

                                for(uint32_t i=0; i<num; i++)
                                {
                                        if (i == largestIdx)
#ifdef NOISY
                                        {
                                                printf("  ~~~~  ");
#endif
                                                continue;
#ifdef NOISY
                                        }
#endif
//                                      if ((a[largestIdx] - a[i]) > 7.0)
                                        if ((a[largestIdx] - a[i]) > 9.0)
                                        {
                                                pulled = true;
                                                a[i] += (double)Global::wave[sNum[i]][syncPt[i]];
#ifdef NOISY
                                                printf(" %6d ", Global::wave[sNum[i]][syncPt[i]]);
#endif
                                                syncPt[i]++;
                                        }
#ifdef NOISY
                                        else
                                                printf("        ");
#endif
                                        // See if we've run out of stream to synthesize...
                                        if (syncPt[i] >= Global::waveLen[sNum[i]])
                                                return swLen;
                                }

#ifdef NOISY
                                printf("\n");

                                for(uint32_t i=0; i<num; i++)
                                        printf("<%6.2lf>", a[i]);
#endif
                                // Calculate the standard deviation at this new point
                                sigma = StdDeviation(a, num);//, &mean);

                                // Calculate the following sigma & lookahead too
                                for(uint32_t i=0; i<num; i++)
                                        b[i] = (double)Global::wave[sNum[i]][syncPt[i]];

                                double nextSigma = StdDeviation(b, num);
                                lookahead = LookaheadWave(syncPt, num);
#ifdef NOISY
                                printf(" ==> %.2lf%s", sigma, (sigma <= SIGMA_THRESHOLD ? " (OK)" : (!pulled ? "(STEP)" : "")));
                                printf(" (lookahead=%.2lf, %d)\n", nextSigma, lookahead);

for(uint32_t j=0; j<num; j++) { printf("{%d}%s", syncPt[j], (j < num - 1 ? "-" : "")); }
printf("\n");//*/
#endif
                                // A handful of heuristics.  If the sigma is OK, we're done
                                if (sigma < 2.5)
                                        keepWorking = false;

                                // If we couldn't pull any streams forward for some reason, punt
                                if (!pulled)
                                        keepWorking = false;

                                // If the average of the current & next sigmas is OK, go ahead
                                if (((sigma + nextSigma) / 2.0) < 2.5)
                                        keepWorking = false;

                                // If the next sigma is OK and the lookahead too, move on
                                if ((nextSigma < 2.5) && (lookahead > 10))
                                        keepWorking = false;
                        }
                        while (keepWorking);

                        // Now, we've re-established sync, now fill the appropriate spots

/*
So the algorithm to do this proppaly is like so:

Step 1: Find the smallest step forward through the section; set amplitude to 0.2.
Step 2: Loop thru the other streams and see if they correspond.  If so, advance their stream start & add 0.2 to the amplitude for each one seen.
Step 3: Add the mean value of the corresponding values to the synthesized waveform + wave amplitude
Step 4: If all of the stream starts are at the end of the section, we are done.
[maybe.  how do you handle the case where you have streams that look like so:

200
150 50
150 50
175 25
200

.......................................|
.............................|.........|
.............................|.........|
..................................|....|
.......................................|

this will fail to add in the correct amount for the times that had no part in the anomaly.  so to fix this, maybe have a baseline value that bumps upward to show where the "cursor" (so to speak) is, so at the end we can subtract the mean that we got above from this baseline to get the remaining time to synthesize, like so:

synthWave[swLen] = (uint32_t)((mean - baseline) + 0.5)
swAmplitude[swLen] = amplitude;

no need for that.  the times that step directly to the end can safely be ignored.
]

150 -> 150 (ignores 175)
150 goes in, ampl. .4, we now have:

200
50
50
175 25
200

only one is left that we can grab:

175

last pulse was at 150 (40%), so we need to subtract that out

-> 25

so maybe we'd do like so:
-1) create a list containing the start values for the start of the anomalous section and replace them as the counters advance.
0) are streams at the end? if so, stuff in the final value to the synthesized wave & stop, otherwise:
1) find the smallest in the current row and use that as a basis.
2) find streams that correlate, if any
3) advance the streams that correlate & subtract out the minimum from the other rows
4) go to 0)

200
150 50
150 50
175 25
200

smallest is 150: pick 150 -> picks up 150
advance the 150 rows & subtract 150 from the rest:

50
50
50
25 25
50

pick smallest: 25 -> picks up nothing
advance the 25 & subtract it from the rest:

25
25
25
25
25

We've reached the end of the section, add the final amount

--------------------------------------------------------------------------------
syncSave before: [2->4][2->6][2->4][2->6][2->6]
syncSave after: [3->4][3->6][3->4][3->6][3->6]

currentRow: [120.89][49.97][46.98][49.98][48.98] => 28.78

2: [3->4][3->6][3->4][3->6][3->6]

a[0]=46.98: [120.89] ==> 36.95
        [49.97]{17.00} ==> 1.50
        <106.00>[49.97][49.98]{28.00} ==> 1.41
        [49.97][49.98][48.98]{26.00} ==> 1.22

3: [3->4][4->6][4->4][4->6][4->6]

a[0]=17.00: [73.91] ==> 28.45
        <56.00>[28.00]{45.00} ==> 5.50
        [28.00][26.00]{47.00} ==> 4.78

4: [3->4][5->6][4->4][5->6][5->6]

a[0]=45.00: [56.91] ==> 5.95
        [56.00]{32.00} ==> 5.50
        <32.00>[56.00][47.00]{32.00} ==> 4.78

5: [3->4][6->6][4->4][6->6][6->6]

a[0]=11.91: <32.00>

6: [4->4][6->6][4->4][6->6][6->6]

 2->4    2->6    2->4    2->6    2->6
[120.89][ 49.97][ 46.98][ 49.98][ 48.98] => 28.78
( 73.91)  17.00  106.00   28.00   26.00
( 56.91)  56.00           45.00   47.00
  32.00   32.00           32.00   32.00

*/
#ifdef NOISY
/*printf("syncSave before: ");
for(uint32_t i=0; i<num; i++)
        printf("[%d->%d]", syncSave[i], syncPt[i]);
printf("\n");//*/

for(uint32_t i=0; i<num; i++)
{
        printf("Sync %d: %d -> %d (", sNum[i], syncSave[i], syncPt[i]);
        for(uint32_t j=syncSave[i]; j<=syncPt[i]; j++)
                printf("%d%s", Global::wave[sNum[i]][j], (j == syncPt[i] ? "" : ", "));
        printf(")\n");
}
#endif
                        // Collect the current row
                        for(uint32_t i=0; i<num; i++)
                        {
                                currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
                                syncSave[i]++;
                        }

/*printf("syncSave after: ");
for(uint32_t i=0; i<num; i++)
        printf("[%d->%d]", syncSave[i], syncPt[i]);
printf("\n");
printf("currentRow: ");
for(uint32_t i=0; i<num; i++)
        printf("[%.2lf]", currentRow[i]);
printf(" => %.2lf\n", StdDeviation(currentRow, num));//*/

//printf("%d: ", swLen);

                        float ampStep = 1.0f / (float)num;

                        while (true)
                        {
                                if (Equal(syncSave, syncPt, num))
                                        break;

                                float amplitude = ampStep;
                                uint32_t smallestIdx = SmallestIndex(currentRow, num);
                                uint32_t cnt = 1;
                                a[0] = currentRow[smallestIdx];
//printf("a[0]=%.2lf: ", a[0]);

                                for(uint32_t i=0; i<num; i++)
                                {
                                        // Make sure we haven't gone beyond the end of any of the streams we're looking at
                                        if (syncSave[i] >= Global::streamLen[i])
                                                return swLen;

                                        if (i == smallestIdx)
                                        {
                                                // We always pick up the smallest in the current set, since that's our basis
                                                currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
                                                syncSave[i]++;
//printf("<%.2lf>", currentRow[i]);
                                                continue;
                                        }

                                        // Save the current row & subtract out the smallest as a basis
                                        double crSave = currentRow[i];
                                        currentRow[i] -= a[0];

                                        // Don't bother with stuff that's out of range (@ current sync point!)  :-P
                                        if (syncSave[i] == syncPt[i])
                                                continue;

                                        // Try to find correlates
                                        a[cnt] = crSave;
                                        cnt++;
                                        sigma = StdDeviation(a, cnt);
//for(uint32_t j=1; j<cnt; j++)
//      printf("[%.2lf]", a[j]);

//relaxed sigma...
                                        if (sigma <= 5.5)//SIGMA_THRESHOLD)
                                        {
                                                // If the sigma is in range for this item, pick it up
                                                currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
                                                syncSave[i]++;

//printf("{%.2lf}", currentRow[i]);
                                                amplitude += ampStep;
                                        }
                                        else
                                        {
                                                // Otherwise reject the last item; it's too far out
                                                cnt--;
                                        }
//printf(" ==> %.2lf\n\t", sigma);
                                }

                                // Get the average of the items picked up and add it to the wave
                                sigma = StdDeviation(a, cnt, &mean);
                                synthWave[swLen] = (uint32_t)(mean + 0.5);
                                swAmplitude[swLen] = amplitude;
                                swLen++;
                        }

                        sigma = StdDeviation(currentRow, num, &mean);
                        synthWave[swLen] = mean;
                        swAmplitude[swLen] = 1.0f;
                        swLen++;
                }
        }

        return swLen;
}
#endif


bool ResyncWave(uint32_t * sync, uint32_t num, int8_t dir = 1);
bool ResyncWave(uint32_t * sync, uint32_t num, int8_t dir/*= 1*/)
{
        double a[6];

        for(uint32_t i=0; i<num; i++)
        {
                a[i] = (double)Global::wave[sNum[i]][sync[i]];
                sync[i] += dir;
        }

#ifdef NOISY
        for(uint32_t i=0; i<num; i++)
                printf("<%4d>", (uint32_t)a[i]);

        printf(" {la: %d}\n", LookaheadWave(sync, num, dir));
#endif

        uint32_t syncCount = 0;
        bool keepGoing = true;

        do
        {
                if (syncCount > 42)
                {
                        printf("Could not synchronize streams...\n");
                        return false;
                }

                for(uint32_t i=0; i<num; i++)
                {
                        if (((dir == 1) && (sync[i] >= Global::waveLen[sNum[i]]))
                                || ((dir == -1) && (sync[i] == 0)))
                        {
                                printf("End of stream...\n");
                                return true;
                        }
                }

                // Deal with the problem.  The heuristic is find the largest value in the set, then advance the other streams by one step until they are "close" to the largest.  If the std. deviation is too large, try again; eventually it should be within good limits.
                uint32_t largestIdx = LargestIndex(a, num);
                bool stepped = false;

                for(uint32_t i=0; i<num; i++)
                {
                        if (i == largestIdx)
#ifdef NOISY
                        {
                                printf(" ~~~~ ");
#endif
                                continue;
#ifdef NOISY
                        }
#endif

                        // If the difference between the highest and this is over 7, add the next time interval and let the sync counter, uh, slip :-)
//                      if ((a[largestIdx] - a[i]) > 7.0)
                        if ((a[largestIdx] - a[i]) > 9.0)
                        {
                                stepped = true;
                                double ttn = (double)Global::wave[sNum[i]][sync[i]];
                                a[i] += ttn;
                                sync[i] += dir;
#ifdef NOISY
                                printf(" %4d ", (uint32_t)ttn);
#endif
                        }
#ifdef NOISY
                        else
                                printf("      ");
#endif
                }

                double sigma = StdDeviation(a, num);
                double nextSigma = StdDeviationWave(sync, num);
                uint32_t lookahead = LookaheadWave(sync, num, dir);
#ifdef NOISY
                printf("\n");

                for(uint32_t i=0; i<num; i++)
                        printf("<%4d>", (uint32_t)a[i]);

                printf(" ==> %.2lf ", sigma);
                printf("(next: %.2lf, la: %d)\n", nextSigma, lookahead);
#endif
                syncCount++;

                // A handful of heuristics...  1st: If sigma is good, we're done
                if (sigma < 2.50)
                        keepGoing = false;

                // Couldn't pull, so punt for now
                if (!stepped)
                        keepGoing = false;

                // If the next sigma + current / 2 is OK, move on
                if (((sigma + nextSigma) / 2.0) < 2.50)
                        keepGoing = false;

                // If the next sigma is good, and the lookahead is good, move on
                // [Slight tweak: only do this if the current sigma isn't terrible]
                if ((nextSigma < 2.50) && (lookahead > 10) && (sigma < 20.0))
                        keepGoing = false;
        }
        while (keepGoing);

        return true;
}


bool AttemptToFindStart(uint32_t * sync, uint32_t num)
{
        // Sanity check...
        for(uint32_t i=0; i<num; i++)
        {
                if (sync[i] == 0)
                        return true;
        }

        bool keepGoing = true;

        do
        {
                double sigma = StdDeviationWave(sync, num);

                if (sigma > 2.50)
                {
                        bool result = ResyncWave(sync, num, -1);

                        if (result == false)
                                return false;
                }
                else
                {
                        for(uint32_t i=0; i<num; i++)
                                sync[i]--;
                }

                for(uint32_t i=0; i<num; i++)
                {
                        if (sync[i] == 0)
                                keepGoing = false;
                }
        }
        while (keepGoing);

//#ifdef NOISY
#if 1
        for(uint32_t i=0; i<num; i++)
                printf("Sync for stream %d: %d\n", i, sync[i]);
        printf("Lookahead: %d\n", LookaheadWave(sync, num));
#endif

        return true;
}


bool InitialSync(uint32_t * sync, uint32_t num)
{
        // Let's try a different approach to this...
        bool lookAt[6] = { true, true, false, false, false, false };

//      while (FindSyncBetweenFirstTwo(sync) == false)
        while (FindSyncBetweenFirstTwo(sync) < 1000)
        {
                for(uint32_t i=0; i<num; i++)
                {
                        sync[i] += 200;

                        if (sync[i] >= Global::waveLen[sNum[i]])
                        {
                                printf("Could not find initial sync!!\n");
                                return false;
                        }
                }
        }

        uint32_t lookahead = LookaheadWave(sync, 2);

        if (lookAt[0])
                printf("Found sync for streams 1 & 2 at %d and %d (la=%d)\n", sync[0], sync[1], lookahead);

        // Add in the others, one at a time and try to synchronize them
        uint32_t syncBest[6];

        for(uint32_t i=2; i<num; i++)
        {
                lookAt[i] = true;
                uint32_t bestLA = 0;

                for(uint32_t j=0; j<200; j++)
                {
                        uint32_t syncSave = sync[i];

                        for(uint32_t k=0; k<200; k++)
                        {
                                uint32_t la = LookaheadWave(sync, i + 1);

                                if (la > bestLA)
                                {
                                        bestLA = la;
                                        memcpy(syncBest, sync, sizeof(syncBest));
                                }

                                sync[i]++;
                        }

                        sync[i] = syncSave;

                        // Kick all the other streams forward to check the rest...
                        for(uint32_t k=0; k<i; k++)
                                sync[k]++;
                }

                memcpy(sync, syncBest, sizeof(syncBest));
        }

        for(uint32_t i=0; i<num; i++)
                printf("Sync for stream %d: %d (%s)\n", i, sync[i], (lookAt[i] ? "yes" : "no"));
        printf("Lookahead: %d\n", LookaheadWave(sync, num));

        return true;
}


uint32_t LoopLookahead(uint32_t trackNum, uint32_t start, uint32_t loopPoint, bool * test = NULL);
uint32_t LoopLookahead(uint32_t trackNum, uint32_t start, uint32_t loopPoint, bool * test/*= NULL*/)
{
        uint32_t count = 0;
        bool equal = true;
        double a[2];

        for(uint32_t i=loopPoint, pos=start; i<Global::synWaveLen[trackNum]; i++, pos++)
        {
                a[0] = Global::synWave[trackNum][pos];
                a[1] = Global::synWave[trackNum][i];
                double sigma = StdDeviation(a, 2);

                if (sigma > 2.50)
                {
                        equal = false;
                        break;
                }

                count++;
        }

        if (test)
                *test = equal;

        return count;
}


#define NOISYSYNTH
//#define DEBUGSYNTH
void SynthesizeTrack(uint32_t trackNum)
{
        uint32_t num = 0;
        uint32_t sync[6] = { 0 };

        for(uint32_t i=0; i<Global::numStreams; i++)
        {
                // Skip streams that aren't on the current track & BITS captures
                if ((Global::stream[i]->location != trackNum)
                        || (Global::stream[i]->captureType == 2))
                        continue;

#ifdef DEBUGSYNTH
printf("sNum[%d] = %d\n", num, i);
#endif
                // Save the stream # so we can find it later
                sNum[num++] = i;
        }

        // We should just bypass the synthesis part of this, not the loop point finding...  :-P  !!! FIX !!!
        if (num == 1)
        {
                printf("%.2f: Single stream, length: %d\n", (float) trackNum / 4.0f, Uint32LE(Global::stream[sNum[0]]->dataLength));

                memcpy(Global::synWave[trackNum], Global::wave[sNum[0]], Global::waveLen[sNum[0]] * 4);
                Global::synWaveLen[trackNum] = Global::waveLen[sNum[0]];
        }

/*
This is failing for track 9.00: it finds an initial sync @ 202, 200, la=1257 but real sync point is 0,2,2 (it says 5,0,0, which is wrong)
*/
        if (num > 1)
        {
                InitialSync(sync, num);
                AttemptToFindStart(sync, num);

                uint32_t count = 0, bloops = 0;
                bool keepGoing = true;
                Global::synWaveLen[trackNum] = 0;
                double mean;

                do
                {
                        double sigma = StdDeviationWave(sync, num, &mean);

                        if (sigma <= 2.50)
                        {
                                count++;
                                Global::synWave[trackNum][Global::synWaveLen[trackNum]] = (uint32_t)(mean + 0.5);
                                Global::synWaveLen[trackNum]++;

                                for(uint32_t i=0; i<num; i++)
                                        sync[i]++;
                        }
                        else
                        {
#ifdef NOISYSYNTH
                        uint32_t startOfSynth = Global::synWaveLen[trackNum];
#endif
#ifdef NOISYSYNTH
                                if (bloops < 2)
                                {
                                        printf("Bloop seen at pos ");
                                        for(uint32_t k=0; k<num; k++)
                                                printf("%d%s", sync[k], (k == num - 1 ? "" : ", "));
                                        printf("...\n");
                                }
#endif
                                bloops++;       // Count the # of times we've had to punt
                                uint32_t syncSave[6] = { 0 };
                                memcpy(syncSave, sync, sizeof(syncSave));
                                bool result = ResyncWave(sync, num);

                                // Now go through the catty wampus crap and zero out the parts that aren't all ones across the board
                                double timeToOnes = 0;
                                double a[6], currentRow[6];

                                // Collect the current row
                                for(uint32_t i=0; i<num; i++)
                                {
                                        currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
                                        syncSave[i]++;
                                }
#ifdef DEBUGSYNTH
for(uint32_t i=0; i<num; i++)
        printf("syncSave[%d]=%d (sync: %d)\n", i, syncSave[i], sync[i]);
#endif

                                while (true)
                                {
                                        if (Equal(syncSave, sync, num))
                                                break;
#ifdef DEBUGSYNTH
printf("syncSave: ");
for(uint32_t i=0; i<num; i++)
        printf("%d%s", syncSave[i], (i == num - 1 ? "" : ", "));
printf("\n    sync: ");
for(uint32_t i=0; i<num; i++)
        printf("%d%s", sync[i], (i == num - 1 ? "" : ", "));
printf("\n");
#endif

                                        uint32_t amplitude = 1;
                                        uint32_t smallestIdx = SmallestIndex(currentRow, num);
                                        uint32_t cnt = 1;
                                        a[0] = currentRow[smallestIdx];
#ifdef DEBUGSYNTH
printf("a[0]=%.2lf:\n", a[0]);
#endif

                                        for(uint32_t i=0; i<num; i++)
                                        {
                                                // Make sure we haven't gone beyond the end of any of the streams we're looking at
                                                if (syncSave[i] >= Global::waveLen[sNum[i]])
                                                {
                                                        printf("Ran past end of stream %d (stream %d): syncSave=%d, waveLen=%d\n", i, sNum[i], syncSave[i], Global::waveLen[sNum[i]]);
                                                        return;
                                                }

                                                if (i == smallestIdx)
                                                {
                                                        // Maybe we should chuck this above the if()?
//                                                      if (syncSave[i] == sync[i])
//                                                              continue;

                                                        // We always pick up the smallest in the current set, since that's our basis
//the reason this works is because this simulates the "add next and subtract smallest", since the subtraction of the smallest would always be zero.  :-P
                                                        currentRow[i] = Global::wave[sNum[i]][syncSave[i]];
//                                                      currentRow[i] += Global::wave[sNum[i]][syncSave[i]];
                                                        syncSave[i]++;
#ifdef DEBUGSYNTH
printf("\t<%.2lf> (i=%d)\n", currentRow[i], i);
#endif
                                                        continue;
                                                }

                                                // Save the current row & subtract out the smallest as a basis
                                                double crSave = currentRow[i];
                                                currentRow[i] -= a[0];

                                                // Don't bother with stuff that's out of range (@ current sync point!)  :-P
                                                if (syncSave[i] == sync[i])
                                                        continue;

                                                // Try to find correlates
                                                a[cnt] = crSave;
                                                cnt++;
                                                sigma = StdDeviation(a, cnt);
#ifdef DEBUGSYNTH
printf("\t");
for(uint32_t j=1; j<cnt; j++)
        printf("[%.2lf]", a[j]);
#endif

//relaxed sigma...
                                                if (sigma <= 5.5)//SIGMA_THRESHOLD)
                                                {
                                                        // If the sigma is in range for this item, pick it up
//                                                      currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
                                                        currentRow[i] += (double)Global::wave[sNum[i]][syncSave[i]];
                                                        syncSave[i]++;

#ifdef DEBUGSYNTH
printf("{%.2lf}", currentRow[i]);
#endif
                                                        amplitude++;
                                                }
                                                else
                                                {
                                                        // Otherwise reject the last item; it's too far out
                                                        cnt--;
                                                }
#ifdef DEBUGSYNTH
printf(" ==> %.2lf\n", sigma);
#endif
                                        }

                                        // Get the average of the items picked up and add it to the wave
                                        sigma = StdDeviation(a, cnt, &mean);
                                        timeToOnes += mean;

                                        // Did we find the line of ones yet?  If so, chuck it in the wave
                                        if (amplitude == num)
                                        {
                                                Global::synWave[trackNum][Global::synWaveLen[trackNum]] = (uint32_t)(timeToOnes + 0.5);
                                                Global::synWaveLen[trackNum]++;
                                                timeToOnes = 0;
                                        }
#ifdef DEBUGSYNTH
printf("\t");
for(uint32_t j=0; j<num; j++)
        printf("_%.2lf_", currentRow[j]);
printf("\n\n");
#endif
                                }

                                // Finally, chuck the final line of ones in
                                sigma = StdDeviation(currentRow, num, &mean);
                                timeToOnes += mean;
                                Global::synWave[trackNum][Global::synWaveLen[trackNum]] = (uint32_t)(timeToOnes + 0.5);
                                Global::synWaveLen[trackNum]++;
#ifdef NOISYSYNTH
printf("Blip @ %d: ", startOfSynth);
for(uint32_t w=startOfSynth; w<Global::synWaveLen[trackNum]; w++)
        printf("%d%s", Global::synWave[trackNum][w], (w == Global::synWaveLen[trackNum] - 1 ? "" : ", "));
printf("\n");
#endif
                        }

                        for(uint32_t i=0; i<num; i++)
                        {
                                if (sync[i] >= Global::waveLen[sNum[i]])
                                {
                                        keepGoing = false;
                                        break;
                                }
                        }
                }
                while (keepGoing);

                printf("\n%.2f: Synthesis done. Bit count: %d, bloop count: %d\n", (float)trackNum / 4.0f, count, bloops);
        }

        // Now, find the loop point and make BITS from it
/*
Now, we take the track and make an ouroborus out of it by finding the sync point.  We can also possibly use it to get rid of more false ones like we did above with the resynchronization code.  So how do we do this?

We could use an approach similar to the original sync between 1 & 2, but there's no guarantee that the lookahead will be good right away.  So we need to do a combination where we count the places where it coincides well and plow through spots where it doesn't (and possibly mask any ones in those areas if it's a good match).

We can start with the estimated loop point that the Applesauce hardware found, by walking it a bit before that point we can make sure we don't miss it.
*/
        // We use the loop point from stream #1 because we're lazy
        uint32_t estLoopTime = Uint32LE(Global::stream[sNum[0]]->estLoopPoint);
        uint32_t curTime = 0;
        uint32_t pos = 0, startOfLoop = 0;
        double a[2];
        uint32_t maxLookahead = 0;
        uint32_t maxPos = 0;
        bool full, full2, full0, full1;
        uint32_t blips, worstBlips, longestStretch, stretch, maxLongestStretch = 0, worstBlipsPos = 0, mlsPos = 0;
        worstBlips = 1e9;

        // Get the index to the correct time in "pos"
        while (curTime < estLoopTime)
                curTime += Global::synWave[trackNum][pos++];

printf("%.2f: estLoopTime=%d, curTime=%d, pos=%d\n", (float)trackNum / 4.0f, estLoopTime, curTime, pos);
        uint32_t posSave = pos;

        // "pos" now holds the start of the search for our loop point.
        for(uint32_t i=0; i<200; i++)
        {
                uint32_t lookahead = 0;
                uint32_t j, pos0;
                bool first = true;
                blips = 0;
                longestStretch = stretch = 0;
                startOfLoop = 0;//shouldn't be necessary...

                for(j=(pos+i-100), pos0=0; j<Global::synWaveLen[trackNum]; j++, pos0++)
                {
                        uint32_t lla0 = LoopLookahead(trackNum, pos0 + 0, j, &full0);
                        uint32_t lla1 = LoopLookahead(trackNum, pos0 + 1, j, &full1);

                        lookahead += lla0;
                        stretch = lla0;

                        if (stretch > longestStretch)
                                longestStretch = stretch;

                        if (full0)
                        {
                                break;
                        }
                        else if (full1 && first) // only check 1st time (for now)
                        {
                                startOfLoop = 1;
                                lookahead = lla1;
                                stretch = lla1;
                                break;
                        }

                        if (first)
                                first = false;

                        pos0 += lla0;
                        j += lla0;

                        a[0] = Global::synWave[trackNum][pos0++];
                        a[1] = Global::synWave[trackNum][j++];
                        double sigma = StdDeviation(a, 2);

#if 0
if (first)
{
        first = false;
//      printf("First blip: %.0lf, %.0lf (sigma=%.2lf) [pos=%d, la=%d]\n", a[0], a[1], sigma, i + pos - 100, pos0);

        if (trackNum == 137)
        {
                uint32_t lla = LoopLookahead(trackNum, 0, j, &full);
                uint32_t llap1 = LoopLookahead(trackNum, 1, j, &full2);
                printf("%.2f: pos0=%d, j=%d, sigma=%lf, la=%d%s, la+1=%d%s", (float)trackNum / 4.0f, pos0, j, sigma, lla, (full ? " *equal" : ""), llap1, (full2 ? " *equal" : ""));
        }
}
#endif

/*                      if (sigma <= 2.50)
                        {
                                lookahead++;
                                stretch++;

                                if (stretch > longestStretch)
                                        longestStretch = stretch;
                        }
                        else//*/
//                      {
                                // Hmm, not in sync.  See if we can re-sync
                                stretch = 0;
                                blips++;

                                do
                                {
                                        if (a[0] < a[1])
                                                a[0] += Global::synWave[trackNum][pos0++];
                                        else
                                                a[1] += Global::synWave[trackNum][j++];

                                        // Maybe put all the other heuristics in too?
                                        if (StdDeviation(a, 2) <= 2.50)
                                                break;

                                        // Now you've gone *too* far...
                                        if (j >= Global::synWaveLen[trackNum])
                                                break;
                                }
                                while (true);
//                      }
                }

#if 0
if (trackNum == 137)
        printf(" blips=%d, la=%d\n", blips, lookahead);
#endif

//maybe we should keep an array of the analyses so we can decide afterward which one is best...
                if (blips < worstBlips)
                {
                        worstBlips = blips;
                        worstBlipsPos = pos + i - 100;
                }

                if (lookahead > maxLookahead)
                {
                        maxLookahead = lookahead;
                        maxPos = pos + i - 100;
                }

                if (longestStretch > maxLongestStretch)
                {
                        maxLongestStretch = longestStretch;
                        mlsPos = pos + i - 100;
                }

                if (blips == 0)
                        break;
        }

#if 0
        bitLen[trackNum] = RenderBitstream(bits[trackNum], &byteLen[trackNum], Global::synWave[trackNum], startOfLoop, maxPos);
#endif

        printf("%.2f: Synthesized wave length: %d\n", (float)trackNum / 4.0f    , Global::synWaveLen[trackNum]);
        printf("%s loop point found at position %d (elp%+d, la=%d, startPos=%d, blips=%d, blipsPos=%d, longest stretch=%d, mlsPos=%d, sol=%d)\n----------------------------------------------------------------------\n", (blips == 0 ? "Perfect" : "Possible"), maxPos, maxPos - posSave, maxLookahead, pos, worstBlips, worstBlipsPos, maxLongestStretch, mlsPos, startOfLoop);
/*
Track 1.00:
36110 is the loop point.  But it fails to find it for some reason. [Because it jumped on the first one it found instead of the longest...]
*/

/*      for(uint32_t i=0; i<maxLookahead+5; i++)
        {
                a[0] = Global::synWave[trackNum][i];
                a[1] = Global::synWave[trackNum][i + maxPos];
                double sigma = StdDeviation(a, 2);
                printf("[%.2lf]", sigma);
        }

        printf("\n");//*/
}


/*
Using a histogram, we can get an idea of the drift from the ones:

                @@
                @@
                @@
                @@
                @@
                @@
                @@
                @@        @@
                @@        @@
                @@        @@
                @@        @@
                @@        @@
                @@   @@   @@
           @@   @@   @@   @@
           @@   @@   @@   @@
      @@   @@   @@   @@   @@
-------------------------------------------------
0000|0584|1792|9611|2627|5382|0482|0058|0013|0000
 28   29   30   31   32   33   34   35   36   37

Peak of the ones occurs at 31, with a secondary at 33. Twos drift:

                @@
                @@        @@
                @@        @@
                @@        @@
                @@        @@
                @@        @@
                @@        @@
                @@   @@   @@
      @@        @@   @@   @@
-------------------------------------------------
0000|0545|0277|5482|1223|5168|0252|0082|0102|0000
 57   58   59   60   61   62   63   64   65   66

This looks like stretched out ones (they might all be spurious...):

      @@
      @@
      @@             @@
      @@             @@
      @@   @@        @@
      @@   @@        @@
      @@   @@        @@
      @@   @@   @@   @@        @@
---------------------------------------
0000|0015|0007|0002|0011|0001|0002|0000
 43   44   45   46   47   48   49   50


*/