From 0e78f44265d96a04b4da949480893ab592469051 Mon Sep 17 00:00:00 2001 From: cnlohr Date: Fri, 20 Jan 2023 02:14:22 -0500 Subject: [PATCH] Add some mega lospec tests. --- attic/Makefile | 14 + .../test_using_square_wave_octave_approach.c | 301 ++++++++++++++ ...sing_square_wave_octave_approach_stipple.c | 328 ++++++++++++++++ attic/test_using_work_selection_heap.c | 367 ++++++++++++++++++ attic/test_using_work_selection_table.c | 321 +++++++++++++++ 5 files changed, 1331 insertions(+) create mode 100644 attic/Makefile create mode 100644 attic/test_using_square_wave_octave_approach.c create mode 100644 attic/test_using_square_wave_octave_approach_stipple.c create mode 100644 attic/test_using_work_selection_heap.c create mode 100644 attic/test_using_work_selection_table.c diff --git a/attic/Makefile b/attic/Makefile new file mode 100644 index 0000000..14a3ad6 --- /dev/null +++ b/attic/Makefile @@ -0,0 +1,14 @@ +PROJS:=test_using_square_wave_octave_approach test_using_work_selection_heap test_using_work_selection_table test_using_square_wave_octave_approach_stipple + +all : $(PROJS) + +CFLAGS:=-I../colorchord2/rawdraw -g +LDFLAGS:=-lGL -lm -lpthread -lX11 + +$(PROJS): %: %.c + gcc -o $@ $^ $(CFLAGS) $(LDFLAGS) + +clean : + rm -rf $(PROJS) + + diff --git a/attic/test_using_square_wave_octave_approach.c b/attic/test_using_square_wave_octave_approach.c new file mode 100644 index 0000000..f632bda --- /dev/null +++ b/attic/test_using_square_wave_octave_approach.c @@ -0,0 +1,301 @@ +/* + An experiment in very, very low-spec ColorChord. This technique foregoes + multiplies. + + Approach 1B --- This variant locks the number of updates for any given + integration to exactly the quadrature size. It was to test theory 1. + + Approach 1 is square_wave_octave_approac_stipple. +*/ + +#include +#include + +#define CNFG_IMPLEMENTATION +#define CNFGOGL +#include "rawdraw_sf.h" + +#include "os_generic.h" + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ); + +int lastx = 200; +int lasty = 1000; + +int main() +{ + #define FSPS 16000 + #define OCTAVES 6 + #define BPERO 24 + #define BASE_FREQ 22.5 + #define QUADRATURE_STEP_DENOMINATOR 16384 + + // Careful: It makes a lot of sense to explore these relationships. + #define SAMPLE_Q 4 + #define MAG_IIR 0 + #define RUNNING_IIR 31 + #define COMPLEX_IIR 2 + + #define TEST_SAMPLES 256 + + + int16_t samples[TEST_SAMPLES]; + int i; + + + CNFGSetup( "Example App", 1024, 768 ); + + // Precomputed Tables + int8_t whichoctave[2< 2048, then perform a quadrature step. + int32_t flipdistance[BPERO]; + int binstothis[BPERO*OCTAVES] = { 0 }; + int nextloopbins[BPERO*OCTAVES] = { 0 }; + for( i = 0; i < BPERO; i++ ) + { + double freq = pow( 2, (double)i / (double)BPERO ) * (BASE_FREQ/2.0); + double pfreq = pow( 2, OCTAVES ) * freq; + double spacing = (FSPS / 2) / pfreq / 4; + + flipdistance[i] = QUADRATURE_STEP_DENOMINATOR * spacing; + // Spacing = "quadrature every X samples" + //printf( "%f %d\n", spacing, flipdistance[i] ); + //flipdistance[i] = QUADRATURE_STEP_DENOMINATOR * (int)(spacing+0.5); + binstothis[i] = (int)(spacing+0.5); + nextloopbins[i] = binstothis[i]; + } + + // This is for timing. Not accumulated data. + int32_t quadrature_timing_last[BPERO*OCTAVES] = { 0 }; + + + uint8_t quadrature_state[BPERO*OCTAVES] = { 0 }; + + uint32_t last_accumulated_value[BPERO*OCTAVES*2] = { 0 }; + + uint32_t octave_timing[OCTAVES] = { 0 }; + + int32_t real_imaginary_running[BPERO*OCTAVES*2] = { 0 }; + + uint32_t sample_accumulator = 0; + + int32_t qcount[BPERO*OCTAVES] = { 0 }; + int32_t magsum[BPERO*OCTAVES] = { 0 }; + + int frameno = 0; + double dLT = OGGetAbsoluteTime(); + int samplenoIn = 0; + int sampleno = 0; + double ToneOmega = 0; + while( CNFGHandleInput() ) + { + CNFGClearFrame(); + + frameno++; + float freq = + //pow( 2, (frameno%600)/100.0 ) * 25; + pow( 2, (lastx)/100.0 ) * lastx; + //101; + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + samples[i] = lasty/5 + sin( ToneOmega ) * 127;// + (rand()%128)-64; + ToneOmega += 1 / (double)FSPS * (double)freq * 3.14159 * 2.0; + } + char cts[1024]; + sprintf( cts, "%f %d", freq, sampleno ); + CNFGColor( 0xffffffff ); + CNFGPenX = 2; + CNFGPenY = 2; + CNFGDrawText( cts, 2 ); + + while( OGGetAbsoluteTime() < dLT + TEST_SAMPLES / (double)FSPS ); + dLT += TEST_SAMPLES / (double)FSPS; + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + sample_accumulator += samples[i]; + int octave = whichoctave[(sampleno)&((2<>RUNNING_IIR) + delta; + real_imaginary_running[last_q_bin] = running; + + int q = ++qcount[binno]; + if( q == SAMPLE_Q ) // Effective Q factor. + { + qcount[binno] = 0; + int newmagR = real_imaginary_running[(binno * 2)]; + int newmagI = real_imaginary_running[(binno * 2)+1]; + + real_imaginary_running[(binno * 2)] = newmagR - (newmagR>>COMPLEX_IIR); + real_imaginary_running[(binno * 2)+1] = newmagI - (newmagI>>COMPLEX_IIR); + + // Super-cheap, non-multiply, approximate complex vector magnitude calculation. + newmagR = (newmagR<0)?-newmagR:newmagR; + newmagI = (newmagI<0)?-newmagI:newmagI; + int newmag = + //sqrt(newmagR*newmagR + newmagI*newmagI ); + newmagR > newmagI ? newmagR + (newmagI>>1) : newmagI + (newmagR>>1); + + int lastmag = magsum[binno]; + magsum[binno] = lastmag - (lastmag>>MAG_IIR) + newmag; + + + quadrature_timing_last[binno] += flipdistance[b]*4; + int nextsteps = -(ocative_time - quadrature_timing_last[binno])/QUADRATURE_STEP_DENOMINATOR; + binstothis[binno] = nextsteps/4; + } + nextloopbins[binno] = binstothis[binno]; + } + } + } + + int lx, ly; + for( i = 0; i < BPERO*OCTAVES; i++ ) + { + CNFGColor( (EHSVtoHEX( (i * 256 / BPERO)&0xff, 255, 255 ) << 8) | 0xff ); + float real = real_imaginary_running[i*2+0]; + float imag = real_imaginary_running[i*2+1]; + float mag = sqrt( real * real + imag * imag ); + + mag = (float)magsum[i] * pow( 2, i / (double)BPERO ); + int y = 768 - ((int)(mag / FSPS * 10) >> MAG_IIR); + if( i ) CNFGTackSegment( i*4, y, lx*4, ly ); + lx = i; ly= y; + //printf( "%d %d\n", real_imaginary_running[i*2+0], real_imaginary_running[i*2+1] ); + } + + CNFGSwapBuffers(); + } +} + + + + + +void HandleKey( int keycode, int bDown ) { } +void HandleButton( int x, int y, int button, int bDown ) { } +void HandleMotion( int x, int y, int mask ) { lastx = x; lasty = y; } +void HandleDestroy() { } + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ) +{ + #define SIXTH1 43 + #define SIXTH2 85 + #define SIXTH3 128 + #define SIXTH4 171 + #define SIXTH5 213 + + uint16_t or = 0, og = 0, ob = 0; + + hue -= SIXTH1; //Off by 60 degrees. + + //TODO: There are colors that overlap here, consider + //tweaking this to make the best use of the colorspace. + + if( hue < SIXTH1 ) //Ok: Yellow->Red. + { + or = 255; + og = 255 - ((uint16_t)hue * 255) / (SIXTH1); + } + else if( hue < SIXTH2 ) //Ok: Red->Purple + { + or = 255; + ob = (uint16_t)hue*255 / SIXTH1 - 255; + } + else if( hue < SIXTH3 ) //Ok: Purple->Blue + { + ob = 255; + or = ((SIXTH3-hue) * 255) / (SIXTH1); + } + else if( hue < SIXTH4 ) //Ok: Blue->Cyan + { + ob = 255; + og = (hue - SIXTH3)*255 / SIXTH1; + } + else if( hue < SIXTH5 ) //Ok: Cyan->Green. + { + og = 255; + ob = ((SIXTH5-hue)*255) / SIXTH1; + } + else //Green->Yellow + { + og = 255; + or = (hue - SIXTH5) * 255 / SIXTH1; + } + + uint16_t rv = val; + if( rv > 128 ) rv++; + uint16_t rs = sat; + if( rs > 128 ) rs++; + + //or, og, ob range from 0...255 now. + //Need to apply saturation and value. + + or = (or * val)>>8; + og = (og * val)>>8; + ob = (ob * val)>>8; + + //OR..OB == 0..65025 + or = or * rs + 255 * (256-rs); + og = og * rs + 255 * (256-rs); + ob = ob * rs + 255 * (256-rs); +//printf( "__%d %d %d =-> %d\n", or, og, ob, rs ); + + or >>= 8; + og >>= 8; + ob >>= 8; + + return or | (og<<8) | ((uint32_t)ob<<16); +} + diff --git a/attic/test_using_square_wave_octave_approach_stipple.c b/attic/test_using_square_wave_octave_approach_stipple.c new file mode 100644 index 0000000..777f613 --- /dev/null +++ b/attic/test_using_square_wave_octave_approach_stipple.c @@ -0,0 +1,328 @@ +/* + An experiment in very, very low-spec ColorChord. This technique foregoes + multiplies. + + Approach 1: + + This approach uses a table, like several colorchord algorithms to only + process one octave each cycle. It then uses a table to decide how often + to process each bin. It performs that bin at 4x the bin's sampling + frequency, in quadrature. So that it performs the +real, +imag, -real, + -imag operations over each cycle. + + You can observe an overtone at the current bin - 1.5 octaves! This is + expected, since, it's the inverse of what the DFT of a square wave would + be. + + That is a minor drawback, but the **major** drawback is that any DC + offset, OR lower frequencies present when computing higher frequencies will + induce a significant ground flutter, and make results really inaccurate. + + This will need to be addressed before the algorithm is ready for prime-time. + + NOTE: To explore: + 1) Consider SAMPLE_Q to possibly use all 4 cycles - though this will + add latency, it will be more "accurate" --> YES! This helps! + 2) Use the cursor to move left and right to sweep tone and up and down + to change DC-bias. This exhibits this algo's shortcoming. + TODO: Can we somehow zero out the DC offset? + + + DISCOVERY: This approach (octave-step in octave) is very sensitive FSPS + + Theories: + * If we force all 4 quadrature value to be the same number of + integration cycles, does that solve it? --> NO! IT GETS MESSY (see Approach 1B) +*/ + +#include +#include + +#define CNFG_IMPLEMENTATION +#define CNFGOGL +#include "rawdraw_sf.h" + +#include "os_generic.h" + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ); + +int lastx = 200; +int lasty = 1000; + +int main() +{ + #define FSPS 22100 + #define OCTAVES 6 + #define BPERO 24 + #define BASE_FREQ 22.5 + #define QUADRATURE_STEP_DENOMINATOR 16384 + + // Careful: It makes a lot of sense to explore these relationships. + #define SAMPLE_Q 4 + #define MAG_IIR 0 + #define RUNNING_IIR 31 + #define COMPLEX_IIR 2 + + #define TEST_SAMPLES 256 + + + int16_t samples[TEST_SAMPLES]; + int i; + + + CNFGSetup( "Example App", 1024, 768 ); + + // Precomputed Tables + int8_t whichoctave[2< 2048, then perform a quadrature step. + int32_t flipdistance[BPERO]; + for( i = 0; i < BPERO; i++ ) + { + double freq = pow( 2, (double)i / (double)BPERO ) * (BASE_FREQ/2.0); + double pfreq = pow( 2, OCTAVES ) * freq; + double spacing = (FSPS / 2) / pfreq / 4; + + flipdistance[i] = QUADRATURE_STEP_DENOMINATOR * spacing; + // Spacing = "quadrature every X samples" + //printf( "%f %d\n", spacing, flipdistance[i] ); + } + + // This is for timing. Not accumulated data. + int32_t quadrature_timing_last[BPERO*OCTAVES] = { 0 }; + uint8_t quadrature_state[BPERO*OCTAVES] = { 0 }; + + uint32_t last_accumulated_value[BPERO*OCTAVES*2] = { 0 }; + + uint32_t octave_timing[OCTAVES] = { 0 }; + + int32_t real_imaginary_running[BPERO*OCTAVES*2] = { 0 }; + + uint32_t sample_accumulator = 0; + + int32_t qcount[BPERO*OCTAVES] = { 0 }; + int32_t magsum[BPERO*OCTAVES] = { 0 }; + + int frameno = 0; + double dLT = OGGetAbsoluteTime(); + int samplenoIn = 0; + int sampleno = 0; + double ToneOmega = 0; + int ops; + while( CNFGHandleInput() ) + { + CNFGClearFrame(); + + frameno++; + float freq = + //pow( 2, (frameno%600)/100.0 ) * 25; + pow( 2, (lastx)/100.0 ) * lastx; + //101; + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + samples[i] = lasty/5 + sin( ToneOmega ) * 127;// + (rand()%128)-64; + ToneOmega += 1 / (double)FSPS * (double)freq * 3.14159 * 2.0; + } + char cts[1024]; + sprintf( cts, "%f %d %f", freq, sampleno, ops/(double)sampleno ); + CNFGColor( 0xffffffff ); + CNFGPenX = 2; + CNFGPenY = 2; + CNFGDrawText( cts, 2 ); + + while( OGGetAbsoluteTime() < dLT + TEST_SAMPLES / (double)FSPS ); + dLT += TEST_SAMPLES / (double)FSPS; + + if( 0 ) + { + memset( real_imaginary_running, 0, sizeof( real_imaginary_running ) ); + memset( last_accumulated_value, 0, sizeof( last_accumulated_value ) ); + memset( quadrature_timing_last, 0, sizeof( quadrature_timing_last ) ); + memset( quadrature_state, 0, sizeof( quadrature_state ) ); + memset( octave_timing, 0, sizeof( octave_timing ) ); + sample_accumulator = 0; + } + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + sample_accumulator += samples[i]; + int octave = whichoctave[(sampleno)&((2< 0 ) + { + ops++; + quadrature_timing_last[binno] += flipdistance[b]; + // This code will get appropriately executed every quadrature update. + + int qstate = quadrature_state[binno] = ( quadrature_state[binno] + 1 ) % 4; + + int last_q_bin = (binno * 2) + ( qstate & 1 ); + int delta = sample_accumulator - last_accumulated_value[last_q_bin]; + last_accumulated_value[last_q_bin] = sample_accumulator; + + if( binno == WATCHBIN ) + printf( "Delta: %d\n", delta ); + + // Qstate = + // (0) = +Cos, (1) = +Sin, (2) = -Cos, (3) = -Sin + if( qstate & 2 ) delta *= -1; + + // Update real and imaginary components with delta. + int running = real_imaginary_running[last_q_bin]; + running = running - (running>>RUNNING_IIR) + delta; + real_imaginary_running[last_q_bin] = running; + + int q = ++qcount[binno]; + if( q == SAMPLE_Q ) // Effective Q factor. + { + qcount[binno] = 0; + int newmagR = real_imaginary_running[(binno * 2)]; + int newmagI = real_imaginary_running[(binno * 2)+1]; + + real_imaginary_running[(binno * 2)] = newmagR - (newmagR>>COMPLEX_IIR); + real_imaginary_running[(binno * 2)+1] = newmagI - (newmagI>>COMPLEX_IIR); + + // Super-cheap, non-multiply, approximate complex vector magnitude calculation. + newmagR = (newmagR<0)?-newmagR:newmagR; + newmagI = (newmagI<0)?-newmagI:newmagI; + int newmag = + //sqrt(newmagR*newmagR + newmagI*newmagI ); + newmagR > newmagI ? newmagR + (newmagI>>1) : newmagI + (newmagR>>1); + + int lastmag = magsum[binno]; + magsum[binno] = lastmag - (lastmag>>MAG_IIR) + newmag; + } + } + } + } + + int lx, ly; + for( i = 0; i < BPERO*OCTAVES; i++ ) + { + CNFGColor( (EHSVtoHEX( (i * 256 / BPERO)&0xff, 255, 255 ) << 8) | 0xff ); + float real = real_imaginary_running[i*2+0]; + float imag = real_imaginary_running[i*2+1]; + float mag = sqrt( real * real + imag * imag ); + + mag = (float)magsum[i] * pow( 2, i / (double)BPERO ); + int y = 768 - ((int)(mag / FSPS * 10) >> MAG_IIR); + if( i ) CNFGTackSegment( i*4, y, lx*4, ly ); + lx = i; ly= y; + //printf( "%d %d\n", real_imaginary_running[i*2+0], real_imaginary_running[i*2+1] ); + } + + CNFGSwapBuffers(); + } +} + + + + + +void HandleKey( int keycode, int bDown ) { } +void HandleButton( int x, int y, int button, int bDown ) { } +void HandleMotion( int x, int y, int mask ) { lastx = x; lasty = y; } +void HandleDestroy() { } + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ) +{ + #define SIXTH1 43 + #define SIXTH2 85 + #define SIXTH3 128 + #define SIXTH4 171 + #define SIXTH5 213 + + uint16_t or = 0, og = 0, ob = 0; + + hue -= SIXTH1; //Off by 60 degrees. + + //TODO: There are colors that overlap here, consider + //tweaking this to make the best use of the colorspace. + + if( hue < SIXTH1 ) //Ok: Yellow->Red. + { + or = 255; + og = 255 - ((uint16_t)hue * 255) / (SIXTH1); + } + else if( hue < SIXTH2 ) //Ok: Red->Purple + { + or = 255; + ob = (uint16_t)hue*255 / SIXTH1 - 255; + } + else if( hue < SIXTH3 ) //Ok: Purple->Blue + { + ob = 255; + or = ((SIXTH3-hue) * 255) / (SIXTH1); + } + else if( hue < SIXTH4 ) //Ok: Blue->Cyan + { + ob = 255; + og = (hue - SIXTH3)*255 / SIXTH1; + } + else if( hue < SIXTH5 ) //Ok: Cyan->Green. + { + og = 255; + ob = ((SIXTH5-hue)*255) / SIXTH1; + } + else //Green->Yellow + { + og = 255; + or = (hue - SIXTH5) * 255 / SIXTH1; + } + + uint16_t rv = val; + if( rv > 128 ) rv++; + uint16_t rs = sat; + if( rs > 128 ) rs++; + + //or, og, ob range from 0...255 now. + //Need to apply saturation and value. + + or = (or * val)>>8; + og = (og * val)>>8; + ob = (ob * val)>>8; + + //OR..OB == 0..65025 + or = or * rs + 255 * (256-rs); + og = og * rs + 255 * (256-rs); + ob = ob * rs + 255 * (256-rs); +//printf( "__%d %d %d =-> %d\n", or, og, ob, rs ); + + or >>= 8; + og >>= 8; + ob >>= 8; + + return or | (og<<8) | ((uint32_t)ob<<16); +} + diff --git a/attic/test_using_work_selection_heap.c b/attic/test_using_work_selection_heap.c new file mode 100644 index 0000000..14524e8 --- /dev/null +++ b/attic/test_using_work_selection_heap.c @@ -0,0 +1,367 @@ +/* + An experiment in very, very low-spec ColorChord. This technique foregoes + multiplies. + + Approach 2: + + Similar approach to Approach 1, in that this uses square waves and quarter + wavelength segments to quadrature encode, but instead of using an octave + at a time, it instead creates a heap to work through every sample. + + That way, the error induced by sample stutter is minimized and the square + waves are as accurate as possible. + + WARNING: With this approach, operations can 'bunch up' so that you may + need to clear many, many ops in a single cycle, so it is not at all + appropirate for being run in an interrupt. + + Another benefit: If sample rate is large, no time is spent working on + samples that don't need work. This is better for a sparse set of ops. + + + TODO: Can we do this approach, but with a fixed table to instruct when to + perform every bin? + + GENERAL OBSERVATION FOR ALL VERSIONS: (applicableto all) If we integrate + only bumps for sin/cos, it seems to have different noise properties. + May be beneficial! +*/ + +#include +#include + +#define CNFG_IMPLEMENTATION +#define CNFGOGL +#include "rawdraw_sf.h" + +#include "os_generic.h" + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ); + +int lastx = 200; +int lasty = 1000; + +#define FSPS 12000 +#define OCTAVES 6 +#define BPERO 24 +#define BASE_FREQ 22.5 +#define QUADRATURE_STEP_DENOMINATOR 16384 + +// Careful: It makes a lot of sense to explore these relationships. +#define SAMPLE_Q 4 +#define MAG_IIR 0 +#define COMPLEX_IIR 2 +#define RUNNING_IIR 31 + +#define TEST_SAMPLES 256 + +int32_t next_heap_events[OCTAVES*BPERO*2] = { 0 }; + +int sineshape; +// This will sort the head node back down into the heap, so the heap will +// remain a min-heap. This is done in log(n) time. But, with our data, it +// experimentally only needs to run for 6.47 iterations per call on average +// assuming 24 BPERO and 6 OCTAVES. + +int heapsteps = 0; +int reheaps = 0; +void PercolateHeap( int32_t current_time ) +{ + reheaps++; + int this_index = 0; + int this_val = next_heap_events[0]; + do + { + heapsteps++; + int left = (this_index * 2 + 1); + int right = (this_index * 2 + 2); + + // At end. WARNING: This heap algorithm is only useful if it's an even number of things. + if( right >= OCTAVES*BPERO ) return; + + int leftval = next_heap_events[left*2]; + int rightval = next_heap_events[right*2]; + int diffleft = leftval - this_val; + int diffright = rightval - this_val; + //printf( "RESORT ID %d / INDEX %d / [%d %d] %d %d %d %d\n", next_heap_events[this_index*2+1], this_index, diffleft, diffright, leftval, rightval, left, right ); + if( diffleft > 0 && diffright > 0 ) + { + // The heap is sorted. We're done. + return; + } + + // Otherwise we have to pick an edge to sort on. + if( diffleft <= diffright ) + { + //printf( "LEFT %d / %d\n", left, this_val ); + int swapevent = next_heap_events[left*2+1]; + next_heap_events[left*2+1] = next_heap_events[this_index*2+1]; + next_heap_events[this_index*2+1] = swapevent; + + next_heap_events[this_index*2+0] = leftval; + next_heap_events[left*2+0] = this_val; + + this_index = left; + } + else + { + //printf( "RIGHT %d\n", right ); + + int swapevent = next_heap_events[right*2+1]; + next_heap_events[right*2+1] = next_heap_events[this_index*2+1]; + next_heap_events[this_index*2+1] = swapevent; + + next_heap_events[this_index*2+0] = rightval; + next_heap_events[right*2+0] = this_val; + + this_index = right; + } + } while( 1 ); +} + + +int main() +{ + + + int16_t samples[TEST_SAMPLES]; + int i; + + + CNFGSetup( "Example App", 1024, 768 ); + + // Make a running counter to count up by this amount every cycle. + // If the new number > 2048, then perform a quadrature step. + int32_t flipdistance[BPERO*OCTAVES]; + for( i = 0; i < BPERO*OCTAVES; i++ ) + { + double freq = pow( 2, (double)i / (double)BPERO ) * (BASE_FREQ/2.0); + double pfreq = freq; + double spacing = (FSPS / 2) / pfreq / 4; + + flipdistance[i] = QUADRATURE_STEP_DENOMINATOR * spacing; + // Spacing = "quadrature every X samples" + + next_heap_events[i*2+1] = i; + } + + // This is for timing. Not accumulated data. + uint8_t quadrature_state[BPERO*OCTAVES] = { 0 }; + + uint32_t last_accumulated_value[BPERO*OCTAVES*2] = { 0 }; + + int32_t real_imaginary_running[BPERO*OCTAVES*2] = { 0 }; + + uint32_t sample_accumulator = 0; + int32_t time_accumulator = 0; + + int32_t qcount[BPERO*OCTAVES] = { 0 }; + int32_t magsum[BPERO*OCTAVES] = { 0 }; + + int frameno = 0; + double dLT = OGGetAbsoluteTime(); + + double ToneOmega = 0; + while( CNFGHandleInput() ) + { + CNFGClearFrame(); + + frameno++; + float freq = + //pow( 2, (frameno%600)/100.0 ) * 25; + pow( 2, (lastx)/100.0 ) * lastx; + //101; + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + samples[i] = lasty/5 + sin( ToneOmega ) * 127;// + (rand()%128)-64; + ToneOmega += 1 / (double)FSPS * (double)freq * 3.14159 * 2.0; + } + char cts[1024]; + sprintf( cts, "%f %d %f %f SINESHAPE: %d", freq, time_accumulator, (double)heapsteps / (double)reheaps, (double)reheaps/(time_accumulator/QUADRATURE_STEP_DENOMINATOR), sineshape ); + CNFGColor( 0xffffffff ); + CNFGPenX = 2; + CNFGPenY = 2; + CNFGDrawText( cts, 2 ); + + while( OGGetAbsoluteTime() < dLT + TEST_SAMPLES / (double)FSPS ); + dLT += TEST_SAMPLES / (double)FSPS; + + + #define WATCHBIN -1 + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + sample_accumulator += samples[i]; + + time_accumulator += QUADRATURE_STEP_DENOMINATOR; + + while( (time_accumulator - next_heap_events[0]) > 0 ) + { + // Event has occurred. + int binno = next_heap_events[1]; + //printf( "%d %d\n", binno, next_heap_events[0] ); + next_heap_events[0] += flipdistance[binno]; + PercolateHeap( time_accumulator ); + + //int j; + //for( j = 0; j < OCTAVES*BPERO; j++ ) printf( "[%d %d]", next_heap_events[j*2+0], next_heap_events[j*2+1] ); + //printf( "\n" ); + + int qstate = quadrature_state[binno] = ( quadrature_state[binno] + 1 ) % 4; + + int last_q_bin = (binno * 2) + ( qstate & 1 ); + int delta; + if( !sineshape ) + { + delta = sample_accumulator - last_accumulated_value[last_q_bin]; + last_accumulated_value[last_q_bin] = sample_accumulator; + } + else + { + // TESTING: Sine Shape - this only integrates bumps for sin/cos + // instead of full triangle waves. + // Observation: For higher frequency bins, this seems to help + // a lot. + // side-benefit, this takes less RAM. + // BUT BUT! It messes with lower frequencies, making them uglier. + delta = sample_accumulator - last_accumulated_value[binno]; + last_accumulated_value[binno] = sample_accumulator; + delta *= 1.4; // Just to normalize the results of the test (not for production) + } + if( binno == WATCHBIN ) + printf( "Delta: %d\n", delta ); + + // Qstate = + // (0) = +Cos, (1) = +Sin, (2) = -Cos, (3) = -Sin + if( qstate & 2 ) delta *= -1; + + // Update real and imaginary components with delta. + int running = real_imaginary_running[last_q_bin]; + running = running - (running>>RUNNING_IIR) + delta; + real_imaginary_running[last_q_bin] = running; + + int q = ++qcount[binno]; + if( q == SAMPLE_Q ) // Effective Q factor. + { + qcount[binno] = 0; + int newmagR = real_imaginary_running[(binno * 2)]; + int newmagI = real_imaginary_running[(binno * 2)+1]; + + real_imaginary_running[(binno * 2)] = newmagR - (newmagR>>COMPLEX_IIR); + real_imaginary_running[(binno * 2)+1] = newmagI - (newmagI>>COMPLEX_IIR); + + // Super-cheap, non-multiply, approximate complex vector magnitude calculation. + newmagR = (newmagR<0)?-newmagR:newmagR; + newmagI = (newmagI<0)?-newmagI:newmagI; + int newmag = + //sqrt(newmagR*newmagR + newmagI*newmagI ); + newmagR > newmagI ? newmagR + (newmagI>>1) : newmagI + (newmagR>>1); + + int lastmag = magsum[binno]; + magsum[binno] = lastmag - (lastmag>>MAG_IIR) + newmag; + } + } + } + + int lx, ly; + for( i = 0; i < BPERO*OCTAVES; i++ ) + { + CNFGColor( (EHSVtoHEX( (i * 256 / BPERO)&0xff, 255, 255 ) << 8) | 0xff ); + float real = real_imaginary_running[i*2+0]; + float imag = real_imaginary_running[i*2+1]; + float mag = sqrt( real * real + imag * imag ); + + mag = (float)magsum[i] * pow( 2, i / (double)BPERO ); + int y = 768 - ((int)(mag / FSPS * 10) >> MAG_IIR); + if( i ) CNFGTackSegment( i*4, y, lx*4, ly ); + lx = i; ly= y; + //printf( "%d %d\n", real_imaginary_running[i*2+0], real_imaginary_running[i*2+1] ); + } + + CNFGSwapBuffers(); + } +} + + + + + +void HandleKey( int keycode, int bDown ) { if( keycode == 'a' && bDown ) sineshape = !sineshape; } +void HandleButton( int x, int y, int button, int bDown ) { } +void HandleMotion( int x, int y, int mask ) { lastx = x; lasty = y; } +void HandleDestroy() { } + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ) +{ + #define SIXTH1 43 + #define SIXTH2 85 + #define SIXTH3 128 + #define SIXTH4 171 + #define SIXTH5 213 + + uint16_t or = 0, og = 0, ob = 0; + + hue -= SIXTH1; //Off by 60 degrees. + + //TODO: There are colors that overlap here, consider + //tweaking this to make the best use of the colorspace. + + if( hue < SIXTH1 ) //Ok: Yellow->Red. + { + or = 255; + og = 255 - ((uint16_t)hue * 255) / (SIXTH1); + } + else if( hue < SIXTH2 ) //Ok: Red->Purple + { + or = 255; + ob = (uint16_t)hue*255 / SIXTH1 - 255; + } + else if( hue < SIXTH3 ) //Ok: Purple->Blue + { + ob = 255; + or = ((SIXTH3-hue) * 255) / (SIXTH1); + } + else if( hue < SIXTH4 ) //Ok: Blue->Cyan + { + ob = 255; + og = (hue - SIXTH3)*255 / SIXTH1; + } + else if( hue < SIXTH5 ) //Ok: Cyan->Green. + { + og = 255; + ob = ((SIXTH5-hue)*255) / SIXTH1; + } + else //Green->Yellow + { + og = 255; + or = (hue - SIXTH5) * 255 / SIXTH1; + } + + uint16_t rv = val; + if( rv > 128 ) rv++; + uint16_t rs = sat; + if( rs > 128 ) rs++; + + //or, og, ob range from 0...255 now. + //Need to apply saturation and value. + + or = (or * val)>>8; + og = (og * val)>>8; + ob = (ob * val)>>8; + + //OR..OB == 0..65025 + or = or * rs + 255 * (256-rs); + og = og * rs + 255 * (256-rs); + ob = ob * rs + 255 * (256-rs); +//printf( "__%d %d %d =-> %d\n", or, og, ob, rs ); + + or >>= 8; + og >>= 8; + ob >>= 8; + + return or | (og<<8) | ((uint32_t)ob<<16); +} + diff --git a/attic/test_using_work_selection_table.c b/attic/test_using_work_selection_table.c new file mode 100644 index 0000000..721455c --- /dev/null +++ b/attic/test_using_work_selection_table.c @@ -0,0 +1,321 @@ +/* + An experiment in very, very low-spec ColorChord. This technique foregoes + multiplies. + + Approach 3: Based on Approach 2, but using a work selection table. + + This won't work for the ch32v003, since the minimum practical table for 6 + octaves at 24BPERO with 12kSPS and bottom freq of 22.5 is about 80kB. +*/ + +#include +#include + +#define CNFG_IMPLEMENTATION +#define CNFGOGL +#include "rawdraw_sf.h" + +#include "os_generic.h" + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ); + +int lastx = 200; +int lasty = 1000; + +#define FSPS 16000 +#define OCTAVES 6 +#define BPERO 24 +#define BASE_FREQ 22.5 +#define QUADRATURE_STEP_DENOMINATOR 16384 + +// Careful: It makes a lot of sense to explore these relationships. +#define SAMPLE_Q 4 +#define MAG_IIR 0 +#define COMPLEX_IIR 2 +#define RUNNING_IIR 31 + +#define TEST_SAMPLES 256 + +int sineshape; +// This will sort the head node back down into the heap, so the heap will +// remain a min-heap. This is done in log(n) time. But, with our data, it +// experimentally only needs to run for 6.47 iterations per call on average +// assuming 24 BPERO and 6 OCTAVES. + + +int main() +{ + + + int16_t samples[TEST_SAMPLES]; + int i; + + + CNFGSetup( "Example App", 1024, 768 ); + + // Not size of table (that's usually larger) but # of samples + // to record the work instructions for. + #define WORKLOOP 6144 + + // Make a running counter to count up by this amount every cycle. + // If the new number > 2048, then perform a quadrature step. + double spacings[BPERO*OCTAVES]; + double runningspace[BPERO*OCTAVES]; + for( i = 0; i < BPERO*OCTAVES; i++ ) + { + double freq = pow( 2, (double)i / (double)BPERO ) * (BASE_FREQ/2.0); + double pfreq = freq; + double spacing = (FSPS / 2) / pfreq / 4; + + // make spacings line up to a denominator of workloop, this makes it + // so you don't get a werid jump at the end of the work loop. + double wdt = WORKLOOP / spacing; + printf( "%f %f %f\n", wdt, spacing, WORKLOOP / (double)((int)wdt+0.5) ); + wdt = (int)(wdt+0.5); + spacing = WORKLOOP / wdt; + spacings[i] = spacing; + } + + + + uint8_t worktable[WORKLOOP*BPERO*OCTAVES]; + int worktablelen = 0; + for( i = 0; i < WORKLOOP; i++ ) + { + int j; + for( j = 0; j < BPERO*OCTAVES; j++ ) + { + if( i >= runningspace[j] ) + { + runningspace[j] += spacings[j]; + worktable[worktablelen++] = j; + } + } + worktable[worktablelen++] = 255; + } + + + // This is for timing. Not accumulated data. + uint8_t quadrature_state[BPERO*OCTAVES] = { 0 }; + + uint32_t last_accumulated_value[BPERO*OCTAVES*2] = { 0 }; + + int32_t real_imaginary_running[BPERO*OCTAVES*2] = { 0 }; + + uint32_t sample_accumulator = 0; + int32_t time_accumulator = 0; + + int32_t qcount[BPERO*OCTAVES] = { 0 }; + int32_t magsum[BPERO*OCTAVES] = { 0 }; + + int frameno = 0; + double dLT = OGGetAbsoluteTime(); + + double ToneOmega = 0; + + int worktableplace = 0; + while( CNFGHandleInput() ) + { + CNFGClearFrame(); + + frameno++; + float freq = + //pow( 2, (frameno%600)/100.0 ) * 25; + pow( 2, (lastx)/100.0 ) * lastx; + //101; + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + samples[i] = lasty/5 + sin( ToneOmega ) * 127;// + (rand()%128)-64; + ToneOmega += 1 / (double)FSPS * (double)freq * 3.14159 * 2.0; + } + char cts[1024]; + sprintf( cts, "%f %d SINESHAPE: %d WT %d", freq, time_accumulator, sineshape , worktablelen); + CNFGColor( 0xffffffff ); + CNFGPenX = 2; + CNFGPenY = 2; + CNFGDrawText( cts, 2 ); + + while( OGGetAbsoluteTime() < dLT + TEST_SAMPLES / (double)FSPS ); + dLT += TEST_SAMPLES / (double)FSPS; + + + #define WATCHBIN -1 + + for( i = 0; i < TEST_SAMPLES; i++ ) + { + sample_accumulator += samples[i]; + + time_accumulator += QUADRATURE_STEP_DENOMINATOR; + + while( 1 ) + { + int wtp = worktable[worktableplace]; + worktableplace = worktableplace + 1; + if( worktableplace >= worktablelen ) worktableplace = 0; + if( wtp == 255 ) break; + + // Event has occurred. + int binno = wtp; + + //int j; + //for( j = 0; j < OCTAVES*BPERO; j++ ) printf( "[%d %d]", next_heap_events[j*2+0], next_heap_events[j*2+1] ); + //printf( "\n" ); + + int qstate = quadrature_state[binno] = ( quadrature_state[binno] + 1 ) % 4; + + int last_q_bin = (binno * 2) + ( qstate & 1 ); + int delta; + if( !sineshape ) + { + delta = sample_accumulator - last_accumulated_value[last_q_bin]; + last_accumulated_value[last_q_bin] = sample_accumulator; + } + else + { + // TESTING: Sine Shape - this only integrates bumps for sin/cos + // instead of full triangle waves. + // Observation: For higher frequency bins, this seems to help + // a lot. + // side-benefit, this takes less RAM. + // BUT BUT! It messes with lower frequencies, making them uglier. + delta = sample_accumulator - last_accumulated_value[binno]; + last_accumulated_value[binno] = sample_accumulator; + delta *= 1.4; // Just to normalize the results of the test (not for production) + } + if( binno == WATCHBIN ) + printf( "Delta: %d\n", delta ); + + // Qstate = + // (0) = +Cos, (1) = +Sin, (2) = -Cos, (3) = -Sin + if( qstate & 2 ) delta *= -1; + + // Update real and imaginary components with delta. + int running = real_imaginary_running[last_q_bin]; + running = running - (running>>RUNNING_IIR) + delta; + real_imaginary_running[last_q_bin] = running; + + int q = ++qcount[binno]; + if( q == SAMPLE_Q ) // Effective Q factor. + { + qcount[binno] = 0; + int newmagR = real_imaginary_running[(binno * 2)]; + int newmagI = real_imaginary_running[(binno * 2)+1]; + + real_imaginary_running[(binno * 2)] = newmagR - (newmagR>>COMPLEX_IIR); + real_imaginary_running[(binno * 2)+1] = newmagI - (newmagI>>COMPLEX_IIR); + + // Super-cheap, non-multiply, approximate complex vector magnitude calculation. + newmagR = (newmagR<0)?-newmagR:newmagR; + newmagI = (newmagI<0)?-newmagI:newmagI; + int newmag = + //sqrt(newmagR*newmagR + newmagI*newmagI ); + newmagR > newmagI ? newmagR + (newmagI>>1) : newmagI + (newmagR>>1); + + int lastmag = magsum[binno]; + magsum[binno] = lastmag - (lastmag>>MAG_IIR) + newmag; + } + } + } + + int lx, ly; + for( i = 0; i < BPERO*OCTAVES; i++ ) + { + CNFGColor( (EHSVtoHEX( (i * 256 / BPERO)&0xff, 255, 255 ) << 8) | 0xff ); + float real = real_imaginary_running[i*2+0]; + float imag = real_imaginary_running[i*2+1]; + float mag = sqrt( real * real + imag * imag ); + + mag = (float)magsum[i] * pow( 2, i / (double)BPERO ); + int y = 768 - ((int)(mag / FSPS * 10) >> MAG_IIR); + if( i ) CNFGTackSegment( i*4, y, lx*4, ly ); + lx = i; ly= y; + //printf( "%d %d\n", real_imaginary_running[i*2+0], real_imaginary_running[i*2+1] ); + } + + CNFGSwapBuffers(); + } +} + + + + + +void HandleKey( int keycode, int bDown ) { if( keycode == 'a' && bDown ) sineshape = !sineshape; } +void HandleButton( int x, int y, int button, int bDown ) { } +void HandleMotion( int x, int y, int mask ) { lastx = x; lasty = y; } +void HandleDestroy() { } + + +uint32_t EHSVtoHEX( uint8_t hue, uint8_t sat, uint8_t val ) +{ + #define SIXTH1 43 + #define SIXTH2 85 + #define SIXTH3 128 + #define SIXTH4 171 + #define SIXTH5 213 + + uint16_t or = 0, og = 0, ob = 0; + + hue -= SIXTH1; //Off by 60 degrees. + + //TODO: There are colors that overlap here, consider + //tweaking this to make the best use of the colorspace. + + if( hue < SIXTH1 ) //Ok: Yellow->Red. + { + or = 255; + og = 255 - ((uint16_t)hue * 255) / (SIXTH1); + } + else if( hue < SIXTH2 ) //Ok: Red->Purple + { + or = 255; + ob = (uint16_t)hue*255 / SIXTH1 - 255; + } + else if( hue < SIXTH3 ) //Ok: Purple->Blue + { + ob = 255; + or = ((SIXTH3-hue) * 255) / (SIXTH1); + } + else if( hue < SIXTH4 ) //Ok: Blue->Cyan + { + ob = 255; + og = (hue - SIXTH3)*255 / SIXTH1; + } + else if( hue < SIXTH5 ) //Ok: Cyan->Green. + { + og = 255; + ob = ((SIXTH5-hue)*255) / SIXTH1; + } + else //Green->Yellow + { + og = 255; + or = (hue - SIXTH5) * 255 / SIXTH1; + } + + uint16_t rv = val; + if( rv > 128 ) rv++; + uint16_t rs = sat; + if( rs > 128 ) rs++; + + //or, og, ob range from 0...255 now. + //Need to apply saturation and value. + + or = (or * val)>>8; + og = (og * val)>>8; + ob = (ob * val)>>8; + + //OR..OB == 0..65025 + or = or * rs + 255 * (256-rs); + og = og * rs + 255 * (256-rs); + ob = ob * rs + 255 * (256-rs); +//printf( "__%d %d %d =-> %d\n", or, og, ob, rs ); + + or >>= 8; + og >>= 8; + ob >>= 8; + + return or | (og<<8) | ((uint32_t)ob<<16); +} +