Skip to content

Commit

Permalink
What's better than an unoptimized DFT? An even more unoptimized DFT! …
Browse files Browse the repository at this point in the history 
…But now the bass goes fast.
  • Loading branch information
@CaiB
CaiB committed Dec 17, 2023
1 parent e0e06c0 commit 96648f6
Show file tree
Hide file tree
Showing 2 changed files with 86 additions and 79 deletions.
158 changes: 80 additions & 78 deletions ColorChord.NET/NoteFinder/ShinNoteFinderDFT.cs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
using ColorChord.NET.API;
using System;
using System.Diagnostics;
using System.Diagnostics.CodeAnalysis;
using System.Runtime.Intrinsics;
using System.Text;
Expand All @@ -14,7 +15,8 @@ public static class ShinNoteFinderDFT
private const uint USHORT_RANGE = ushort.MaxValue + 1;

private const ushort SINE_TABLE_90_OFFSET = 8;
private static readonly Vector256<short> SinWave = Vector256.Create(0, 3196, 6270, 9102, 11585, 13622, 15136, 16068, 16383, 16068, 15136, 13622, 11585, 9102, 6270, 3196);
//private static readonly Vector256<short> SinWave = Vector256.Create(0, 3196, 6270, 9102, 11585, 13622, 15136, 16068, 16383, 16068, 15136, 13622, 11585, 9102, 6270, 3196);
private static readonly short[] SinWave = new short[] { 0, 3196, 6270, 9102, 11585, 13622, 15136, 16068, 16383, 16068, 15136, 13622, 11585, 9102, 6270, 3196 };

// The above generated by:
/*
Expand All @@ -39,7 +41,7 @@ public static class ShinNoteFinderDFT
/// <summary> The sample rate of the incoming audio signal, and our reference waveforms. </summary>
public static uint SampleRate = 48000;

public static float StartFrequency = 55f;
public static float StartFrequency = 27.5f;

private static uint GlobalSampleCounter = 0;

Expand All @@ -52,24 +54,20 @@ public static class ShinNoteFinderDFT
private static float StartFrequencyOfTopOctave => StartFrequency * MathF.Pow(2, OctaveCount - 1);

/// <summary> Where raw and resampled audio data is stored. </summary>
/// <remarks> Indexed by [Octave][Sample], size is [<see cref="OctaveCount"/>][<see cref="MaxPresentWindowSize"/>] </remarks>
private static short[][] AudioBuffer;
/// <remarks> Size is [<see cref="MaxPresentWindowSize"/>] </remarks>
private static short[] AudioBuffer;

/// <summary> How large the audio buffer should be treated as being, for each bin. </summary>
/// <remarks> Indexed by [Bin], size is [<see cref="BinCount"/>] </remarks>
private static uint[] AudioBufferSizes;

/// <summary> Where in the audio buffer each bin has added data up to. </summary>
/// <remarks> Indexed by [Octave], size is [<see cref="OctaveCount"/>] </remarks>
private static ushort[] AudioBufferAddHeads;
/// <summary> Up to where in the audio buffer data has been added to. </summary>
private static ushort AudioBufferAddHead;

/// <summary> Where in the audio buffer each bin has removed data up to. </summary>
/// <remarks> Indexed by [Bin], size is [<see cref="BinCount"/>] </remarks>
private static ushort[] AudioBufferSubHeads;

/// <summary> Where raw audio data gets stored as it gets resampled down for lower octaves, size is [<see cref="OctaveCount"/>]</summary>
private static int[] InterBufferAccumulators;

/// <summary> How far forward in the sine table this bin should step with every added sample, such that one full sine wave (wrap back to 0) occurs after the number of steps corresponding to the bin frequency. Format is fixed-point 5b+11b. </summary>
/// <remarks> Indexed by [Bin], size is [<see cref="BinsPerOctave"/>] </remarks>
private static DualU16[] SinTableStepSize;
Expand Down Expand Up @@ -100,71 +98,77 @@ static ShinNoteFinderDFT()
[MemberNotNull(
nameof(AudioBuffer),
nameof(AudioBufferSizes),
nameof(AudioBufferAddHeads), nameof(AudioBufferSubHeads),
nameof(InterBufferAccumulators),
nameof(AudioBufferSubHeads),
nameof(SinTableStepSize), nameof(SinTableLocationAdd), nameof(SinTableLocationSub),
nameof(SinProductAccumulators), nameof(CosProductAccumulators),
nameof(RawBinMagnitudes)
)]
public static void Reconfigure()
{
AudioBuffer = new short[OctaveCount][];
AudioBufferSizes = new uint[BinCount];
AudioBufferAddHeads = new ushort[OctaveCount];
AudioBufferSubHeads = new ushort[BinCount];
InterBufferAccumulators = new int[OctaveCount];
SinTableStepSize = new DualU16[BinsPerOctave];
SinTableStepSize = new DualU16[BinCount];
SinTableLocationAdd = new DualU16[BinCount];
SinTableLocationSub = new DualU16[BinCount];
SinProductAccumulators = new DualI64[BinCount];
CosProductAccumulators = new DualI64[BinCount];
RawBinMagnitudes = new float[BinCount];

float TopStart = StartFrequencyOfTopOctave;
//float TopStart = StartFrequencyOfTopOctave;

uint MaxAudioBufferSize = 0;
// Operations that occur on all bins
for (uint Bin = 0; Bin < BinCount; Bin++)
{
uint WrappedBinIndex = Bin % BinsPerOctave;
//uint WrappedBinIndex = Bin % BinsPerOctave;

float TopOctaveBinFreq = CalculateNoteFrequency(TopStart, BinsPerOctave, WrappedBinIndex);
float TopOctaveNextBinFreq = CalculateNoteFrequency(TopStart, BinsPerOctave, WrappedBinIndex + 2);
float IdealWindowSize = WindowSizeForBinWidth(TopOctaveNextBinFreq - TopOctaveBinFreq); // TODO: Add scale factor to shift this from no overlap to -3dB point
float ThisOctaveStart = StartFrequency * MathF.Pow(2, Bin / BinsPerOctave);
float TopOctaveBinFreq = CalculateNoteFrequency(StartFrequency, BinsPerOctave, Bin);
float TopOctaveNextBinFreq = CalculateNoteFrequency(StartFrequency, BinsPerOctave, Bin + 2);
//float IdealWindowSize = WindowSizeForBinWidth(TopOctaveNextBinFreq - TopOctaveBinFreq); // TODO: Add scale factor to shift this from no overlap to -3dB point
uint ThisBufferSize = RoundedWindowSizeForBinWidth(TopOctaveNextBinFreq - TopOctaveBinFreq, TopOctaveBinFreq, SampleRate);
//ushort ThisBufferSize = (ushort)Math.Ceiling(IdealWindowSize);
AudioBufferSizes[Bin] = Math.Min(MAX_WINDOW_SIZE, ThisBufferSize);
AudioBufferSizes[Bin] = ThisBufferSize; //Math.Min(MAX_WINDOW_SIZE, ThisBufferSize);

if (Bin == 0) { MaxAudioBufferSize = ThisBufferSize; }
MaxAudioBufferSize = Math.Max(MaxAudioBufferSize, ThisBufferSize);

AudioBufferSubHeads[Bin] = (ushort)(WrappedBinIndex == 0 ? 0 : (1 - ThisBufferSize + MaxAudioBufferSize) % MaxAudioBufferSize);


//uint BinInTop = StartOfTopOctave + Bin;
float NCOffset = SampleRate / (AudioBufferSizes[Bin] * 2F);
float StepSizeNCL = USHORT_RANGE * (CalculateNoteFrequency(StartFrequency, BinsPerOctave, Bin) - NCOffset) / SampleRate;
float StepSizeNCR = USHORT_RANGE * (CalculateNoteFrequency(StartFrequency, BinsPerOctave, Bin) + NCOffset) / SampleRate;
SinTableStepSize[Bin].NCLeft = (ushort)Math.Round(StepSizeNCL);
SinTableStepSize[Bin].NCRight = (ushort)Math.Round(StepSizeNCR);

}
MaxPresentWindowSize = MaxAudioBufferSize;

AudioBuffer = new short[MaxAudioBufferSize];

for (uint Bin = 0; Bin < BinCount; Bin++)
{
//AudioBufferSubHeads[Bin] = (ushort)(Bin == 0 ? 0 : (1 - AudioBufferSizes[Bin] + MaxAudioBufferSize) % MaxAudioBufferSize);
AudioBufferSubHeads[Bin] = (ushort)((1 - AudioBufferSizes[Bin] + MaxAudioBufferSize) % MaxAudioBufferSize);
}

// Operations that occur on only one octave's worth of bins
for (uint Bin = 0; Bin < BinsPerOctave; Bin++)
/*for (uint Bin = 0; Bin < BinsPerOctave; Bin++)
{
uint BinInTop = StartOfTopOctave + Bin;
float NCOffset = SampleRate / (AudioBufferSizes[BinInTop] * 2F);
float StepSizeNCL = USHORT_RANGE * (CalculateNoteFrequency(StartFrequencyOfTopOctave, BinsPerOctave, Bin) - NCOffset) / SampleRate;
float StepSizeNCR = USHORT_RANGE * (CalculateNoteFrequency(StartFrequencyOfTopOctave, BinsPerOctave, Bin) + NCOffset) / SampleRate;
SinTableStepSize[Bin].NCLeft = (ushort)Math.Round(StepSizeNCL);
SinTableStepSize[Bin].NCRight = (ushort)Math.Round(StepSizeNCR);
}
}*/

// All bins again, but needs data calculated after the previous one
for (uint Bin = 0; Bin < BinCount; Bin++) // TODO: See if this can be optimized, low priority since it should happen very rarely
{
uint WrappedBinIndex = Bin % BinsPerOctave;
SinTableLocationSub[Bin].NCLeft = (ushort)(-(SinTableStepSize[WrappedBinIndex].NCLeft * (AudioBufferSizes[Bin] - (WrappedBinIndex == 0 ? 0 : 1))));
SinTableLocationSub[Bin].NCRight = (ushort)(-(SinTableStepSize[WrappedBinIndex].NCRight * (AudioBufferSizes[Bin] - (WrappedBinIndex == 0 ? 0 : 1))));
}

// Operations that occur for each octave
for (int Octave = 0; Octave < OctaveCount; Octave++)
{
AudioBuffer[Octave] = new short[MaxAudioBufferSize];
//uint WrappedBinIndex = Bin % BinsPerOctave;
SinTableLocationSub[Bin].NCLeft = (ushort)(-(SinTableStepSize[Bin].NCLeft * (AudioBufferSizes[Bin] - (Bin == 0 ? 0 : 1))));
SinTableLocationSub[Bin].NCRight = (ushort)(-(SinTableStepSize[Bin].NCRight * (AudioBufferSizes[Bin] - (Bin == 0 ? 0 : 1))));
}

Log.Debug(nameof(ShinNoteFinder) + " bin window lengths:");
Expand Down Expand Up @@ -193,7 +197,7 @@ public static void AddAudioData(ReadOnlySpan<short> newData)
// TODO: This handles the top octave only
for (int i = 0; i < newData.Length; i++)
{
AddAudioDataToOctave(newData[i], OctaveCount - 1);
AddAudioDataToOctave(newData[i]);
GlobalSampleCounter++;
}
}
Expand All @@ -203,60 +207,53 @@ public static void AddAudioData(ReadOnlySpan<short> newData)
for (int i = 0; i < newData.Length; i++)
{
short NewData = (short)(newData[i] * short.MaxValue);
AddAudioDataToOctave(NewData, OctaveCount - 1);
AddAudioDataToOctave(NewData);
GlobalSampleCounter++;
}
}

private static void AddAudioDataToOctave(short newData, int octave)
private static void AddAudioDataToOctave(short newData)
{
int OctaveBinOffset = octave * BinsPerOctave;

ushort BinCount = ShinNoteFinderDFT.BinCount;
// Subtract old data from accumulators
for (int Bin = 0; Bin < BinsPerOctave; Bin++)
for (int Bin = 0; Bin < BinCount; Bin++)
{
int FullBinIndex = OctaveBinOffset + Bin;

// Find where we are in the sine table
DualU16 SinTableLoc = SinTableLocationSub[FullBinIndex];
DualU16 SinTableLoc = SinTableLocationSub[Bin];
DualI16 SinValue = GetSine(SinTableLoc, false);
DualI16 CosValue = GetSine(SinTableLoc, true);

// Multiply the outgoing sample by the correct sine sample
short OldBufferData = AudioBuffer[octave][AudioBufferSubHeads[FullBinIndex]];
short OldBufferData = AudioBuffer[AudioBufferSubHeads[Bin]];
int OldSinProductNCL = SinValue.NCLeft * OldBufferData;
int OldSinProductNCR = SinValue.NCRight * OldBufferData;
int OldCosProductNCL = CosValue.NCLeft * OldBufferData;
int OldCosProductNCR = CosValue.NCRight * OldBufferData;

// Remove the product from the accumulators
SinProductAccumulators[FullBinIndex].NCLeft -= OldSinProductNCL;
SinProductAccumulators[FullBinIndex].NCRight -= OldSinProductNCR;
CosProductAccumulators[FullBinIndex].NCLeft -= OldCosProductNCL;
CosProductAccumulators[FullBinIndex].NCRight -= OldCosProductNCR;
SinProductAccumulators[Bin].NCLeft -= OldSinProductNCL;
SinProductAccumulators[Bin].NCRight -= OldSinProductNCR;
CosProductAccumulators[Bin].NCLeft -= OldCosProductNCL;
CosProductAccumulators[Bin].NCRight -= OldCosProductNCR;

// Advance the buffer and sine table locations
AudioBufferSubHeads[FullBinIndex] = (ushort)((AudioBufferSubHeads[FullBinIndex] + 1) % MaxPresentWindowSize);
AudioBufferSubHeads[Bin] = (ushort)((AudioBufferSubHeads[Bin] + 1) % MaxPresentWindowSize);

DualU16 SinTableStep = SinTableStepSize[Bin];
SinTableLocationSub[FullBinIndex].NCLeft += SinTableStep.NCLeft;
SinTableLocationSub[FullBinIndex].NCRight += SinTableStep.NCRight;

//if (Bin == 6) { Console.Write($"{SubHeadBefore},{SinTableLoc.NCRight},{OldSinProductNCL},"); }
SinTableLocationSub[Bin].NCLeft += SinTableStep.NCLeft;
SinTableLocationSub[Bin].NCRight += SinTableStep.NCRight;
}

// Write new data
ushort HeadBefore = AudioBufferAddHeads[octave];
AudioBuffer[octave][HeadBefore] = newData;
AudioBufferAddHeads[octave] = (ushort)((HeadBefore + 1) % MaxPresentWindowSize);
ushort HeadBefore = AudioBufferAddHead;
AudioBuffer[HeadBefore] = newData;
AudioBufferAddHead = (ushort)((HeadBefore + 1) % MaxPresentWindowSize);

// Add new data to accumulators
for (int Bin = 0; Bin < BinsPerOctave; Bin++)
for (int Bin = 0; Bin < BinCount; Bin++)
{
int FullBinIndex = OctaveBinOffset + Bin;

// Find where we are in the sine table
DualU16 SinTableLoc = SinTableLocationAdd[FullBinIndex];
DualU16 SinTableLoc = SinTableLocationAdd[Bin];
DualI16 SinValue = GetSine(SinTableLoc, false);
DualI16 CosValue = GetSine(SinTableLoc, true);

Expand All @@ -267,25 +264,16 @@ private static void AddAudioDataToOctave(short newData, int octave)
int NewCosProductNCR = CosValue.NCRight * newData;

// Add the product to the accumulators
SinProductAccumulators[FullBinIndex].NCLeft += NewSinProductNCL;
SinProductAccumulators[FullBinIndex].NCRight += NewSinProductNCR;
CosProductAccumulators[FullBinIndex].NCLeft += NewCosProductNCL;
CosProductAccumulators[FullBinIndex].NCRight += NewCosProductNCR;
SinProductAccumulators[Bin].NCLeft += NewSinProductNCL;
SinProductAccumulators[Bin].NCRight += NewSinProductNCR;
CosProductAccumulators[Bin].NCLeft += NewCosProductNCL;
CosProductAccumulators[Bin].NCRight += NewCosProductNCR;

// Advance the sine table locations
DualU16 SinTableStep = SinTableStepSize[Bin];
SinTableLocationAdd[FullBinIndex].NCLeft += SinTableStep.NCLeft;
SinTableLocationAdd[FullBinIndex].NCRight += SinTableStep.NCRight;

//if (Bin == 6) { Console.WriteLine($"{HeadBefore},{SinTableLoc.NCRight},{NewSinProductNCL},{SinProductAccumulators[FullBinIndex].NCLeft}"); }
SinTableLocationAdd[Bin].NCLeft += SinTableStep.NCLeft;
SinTableLocationAdd[Bin].NCRight += SinTableStep.NCRight;
}

if (octave > 0 && ((GlobalSampleCounter >> (OctaveCount - 1 - octave)) & 1) != 0) // Need to push audio data down
{
//Console.WriteLine($"At {GlobalSampleCounter}, also processing octave {octave - 1}");
AddAudioDataToOctave((short)((newData + InterBufferAccumulators[octave]) >> 1), octave - 1);
}
else { InterBufferAccumulators[octave] = newData; }
}

public static void CalculateOutput()
Expand Down Expand Up @@ -313,6 +301,7 @@ public static void CalculateOutput()
//double SimpleSq = ((double)SinProductAccumulators[Bin].NCLeft * SinProductAccumulators[Bin].NCLeft) + ((double)CosProductAccumulators[Bin].NCLeft * CosProductAccumulators[Bin].NCLeft);
//RawBinMagnitudes[Bin] = (float)Math.Sqrt(Math.Max(0, SimpleSq)) / AudioBufferSizes[Bin];


float OutBinVal = MathF.Sqrt(RawBinMagnitudes[Bin]) / 8000;
ShinNoteFinder.OctaveBinValues[Bin % BinsPerOctave] += OutBinVal / OctaveCount;
ShinNoteFinder.AllBinValues[Bin] = OutBinVal;
Expand Down Expand Up @@ -341,6 +330,16 @@ public static DualI16 GetSine(DualU16 sineTablePosition, bool shiftForCos) // TO
short ValueUpperL = (short)(((AdjacentLocationL >> 4) == 1 ? -1 : 1) * SinWave[AdjacentLocationL & 0b1111]);
short ValueUpperR = (short)(((AdjacentLocationR >> 4) == 1 ? -1 : 1) * SinWave[AdjacentLocationR & 0b1111]);

//int LocationModifierL = (WholeLocationL << 27) >> 31;
//int LocationModifierR = (WholeLocationR << 27) >> 31;
//short ValueLowerL = (short)(((WholeLocationL ^ LocationModifierL) - LocationModifierL) * SinWave[WholeLocationL & 0b1111]); // The sine vector only contains the positive half of the wave, indices in the range 16...31 are just 0...15 but negative value
//short ValueLowerR = (short)(((WholeLocationR ^ LocationModifierR) - LocationModifierR) * SinWave[WholeLocationR & 0b1111]);

//int AdjLocationModifierL = (AdjacentLocationL << 27) >> 31;
//int AdjLocationModifierR = (AdjacentLocationR << 27) >> 31;
//short ValueUpperL = (short)(((AdjacentLocationL ^ AdjLocationModifierL) - AdjLocationModifierL) * SinWave[AdjacentLocationL & 0b1111]);
//short ValueUpperR = (short)(((AdjacentLocationL ^ AdjLocationModifierR) - AdjLocationModifierR) * SinWave[AdjacentLocationR & 0b1111]);

short FractionalPartL = (short)((((ValueUpperL - ValueLowerL) << 8) * ((sineTablePosition.NCLeft >> 3) & 0xFF)) >> 16); // Wasting the bottom 3 bits of position precision, but otherwise multiplication may overflow int
short FractionalPartR = (short)((((ValueUpperR - ValueLowerR) << 8) * ((sineTablePosition.NCRight >> 3) & 0xFF)) >> 16);
return new() { NCLeft = (short)(ValueLowerL + FractionalPartL), NCRight = (short)(ValueLowerR + FractionalPartR) };
Expand All @@ -365,14 +364,17 @@ public struct DualI64
// Everyone loves magic numbers :)
// These were determined through simulations and regressions, which can be found in the Simulations folder in the root of the ColorChord.NET repository.
private static float BinWidthAtWindowSize(float windowSize) => 50222.5926786413F / (windowSize + 11.483904495504245F);
private static float WindowSizeForBinWidth(float binWidth) => Math.Min(MAX_WINDOW_SIZE, Math.Max(MIN_WINDOW_SIZE, (50222.5926786413F / binWidth) - 11.483904495504245F));
private static float WindowSizeForBinWidth(float binWidth) => (50222.5926786413F / binWidth) - 11.483904495504245F;

private static uint RoundedWindowSizeForBinWidth(float binWidth, float frequency, float sampleRate)
{
float IdealWindowSize = WindowSizeForBinWidth(binWidth);
float PeriodInSamples = sampleRate / frequency;
float PeriodsInWindow = IdealWindowSize / PeriodInSamples;
return (uint)MathF.Round(MathF.Round(PeriodsInWindow) * PeriodInSamples);
float MaxWindowSizePeriods = MathF.Ceiling(MAX_WINDOW_SIZE / PeriodInSamples);
float MinWindowSizePeriods = MathF.Floor(MIN_WINDOW_SIZE / PeriodInSamples);

return (uint)MathF.Round(MathF.Max(MinWindowSizePeriods, MathF.Min(MaxWindowSizePeriods, MathF.Round(PeriodsInWindow))) * PeriodInSamples);
}

private static float CalculateNoteFrequency(float octaveStart, uint binsPerOctave, uint binIndex) => octaveStart * GetNoteFrequencyMultiplier(binsPerOctave, binIndex);
Expand Down
7 changes: 6 additions & 1 deletion ColorChordTests/ShinNoteFinderTest.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,11 @@ namespace ColorChordTests;
[TestClass]
public class ShinNoteFinderTest
{
public ShinNoteFinderTest()
{
ShinNoteFinder a = new("Testing NoteFinder", new());
}

[TestMethod]
public void TestInit()
{
Expand Down Expand Up @@ -109,7 +114,7 @@ public void OutputResponseData()

string[] OutputLines = new string[Steps];

int WindowSize = ShinNoteFinderDFT.MaxPresentWindowSize;
uint WindowSize = ShinNoteFinderDFT.MaxPresentWindowSize;
for (int f = 0; f < Steps; f++)
{
float FreqHere = MinFreq + ((float)f / Steps * (MaxFreq - MinFreq)); // Linear, maybe make log later on?
Expand Down

0 comments on commit 96648f6

Please sign in to comment.