calcBiquad returns nonsense for FILTER_HISHELF case

[jimlenz2000]

Description

in control_sgtl5000.cpp calcBiquad subroutine, the FILTER_HISHELF case needs a "break" statement to prevent falling through to the default code.

Steps To Reproduce Problem

the sketch below returns the same nonsense coefficients regards of HiShelfGain value.

Hardware & Software

Teensy 3.2 and Teensy Audio Library
Audio Workshop Kit
Arduino IDE version: 1.8.10
Teensyduino version: 1.48
Operating system & version: Windows 10

Arduino Sketch

#include <Audio.h>
int coef[5];
void setup() {
  Serial.begin(9600);
  while (!Serial) {} // wait for serial port to connect.
  float HiShelfGain = -6.02;
  calcBiquad(FILTER_HISHELF, 440., HiShelfGain, 0.7071, 524288, 44100, coef);
  for (int i = 0; i < 5; i++) Serial.println(float(coef[i])/524288, 10);
}
void loop() {}

Errors or Incorrect Output

0.25
0.00
0.00
0.00
0.00

[grahamwhaley]

Indeed, that looks broken to me (I happen to be staring at that code, thought I'd confirm it). For ref, snippet from the code:

  case FILTER_LOSHELF:
    b0 = A * ((A+1.0F) - ((A-1.0F)*cosw) + (beta*sinw));
    b1 = 2.0F * A * ((A-1.0F) - ((A+1.0F)*cosw));
    b2 = A * ((A+1.0F) - ((A-1.0F)*cosw) - (beta*sinw));
    a0 = (A+1.0F) + ((A-1.0F)*cosw) + (beta*sinw);
    a1 = 2.0F * ((A-1.0F) + ((A+1.0F)*cosw));
    a2 = -((A+1.0F) + ((A-1.0F)*cosw) - (beta*sinw));
  break;
  case FILTER_HISHELF:
    b0 = A * ((A+1.0F) + ((A-1.0F)*cosw) + (beta*sinw));
    b1 = -2.0F * A * ((A-1.0F) + ((A+1.0F)*cosw));
    b2 = A * ((A+1.0F) + ((A-1.0F)*cosw) - (beta*sinw));
    a0 = (A+1.0F) - ((A-1.0F)*cosw) + (beta*sinw);
    a1 = -2.0F * ((A-1.0F) - ((A+1.0F)*cosw));
    a2 = -((A+1.0F) - ((A-1.0F)*cosw) - (beta*sinw));
  default:
    b0 = 0.5;
    b1 = 0.0;
    b2 = 0.0;
    a0 = 1.0;
    a1 = 0.0;
    a2 = 0.0;

[robsoles]

Hi there. It looks like the calc_biquad I contributed. By memory, been a while, I linked to my source[s] and I tested each case in the switch using a YADDS to sweep the results without becoming upset about what happened - check the source I used, perhaps I erred in my re-arrangement [if any] of their code however I do the feel the urge to point out that these filters have relatively basic mathematical limits and sending it parameters that will return garbage because 'bounds' were broken is hardly impossible. I look forward to seeing the demonstration that it is actually broken, and the fix some genius comes up with - feeling at least a tad responsible, perhaps if you convince me I stuffed it up I will be the genius... Cheers, Rob.

…

On Tue, Dec 29, 2020 at 2:59 AM Graham Whaley ***@***.***> wrote: Indeed, that looks broken to me (I happen to be staring at that code, thought I'd confirm it). For ref, snippet from the code <https://github.com/PaulStoffregen/Audio/blob/master/control_sgtl5000.cpp#L996-L1009> : case FILTER_LOSHELF: b0 = A * ((A+1.0F) - ((A-1.0F)*cosw) + (beta*sinw)); b1 = 2.0F * A * ((A-1.0F) - ((A+1.0F)*cosw)); b2 = A * ((A+1.0F) - ((A-1.0F)*cosw) - (beta*sinw)); a0 = (A+1.0F) + ((A-1.0F)*cosw) + (beta*sinw); a1 = 2.0F * ((A-1.0F) + ((A+1.0F)*cosw)); a2 = -((A+1.0F) + ((A-1.0F)*cosw) - (beta*sinw)); break; case FILTER_HISHELF: b0 = A * ((A+1.0F) + ((A-1.0F)*cosw) + (beta*sinw)); b1 = -2.0F * A * ((A-1.0F) + ((A+1.0F)*cosw)); b2 = A * ((A+1.0F) + ((A-1.0F)*cosw) - (beta*sinw)); a0 = (A+1.0F) - ((A-1.0F)*cosw) + (beta*sinw); a1 = -2.0F * ((A-1.0F) - ((A+1.0F)*cosw)); a2 = -((A+1.0F) - ((A-1.0F)*cosw) - (beta*sinw)); default: b0 = 0.5; b1 = 0.0; b2 = 0.0; a0 = 1.0; a1 = 0.0; a2 = 0.0; — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#343 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABSHQUQUBA6M3G22OGUPKWDSXCTOLANCNFSM4NCTEOOA> .

[grahamwhaley]

Hi @robsoles - thanks for responding. The bug isn't a problem with the calculations themselves, just with a missing break; at the end of the hishelf switch case, which means it falls through and wipes out all the calculated coeffs with the default (benign passthrough) ones... I knocked up a simple example ino to show it (as generally requested in the Issue template):

#include <Audio.h>

void setup() {
  int coeffs[5];

  delay(1000);
  Serial.begin(115200);
  Serial.println("Starting");
  calcBiquad(FILTER_HISHELF, 1000, 0.0, 0.707, 524288, AUDIO_SAMPLE_RATE_EXACT, coeffs);

  for( int i=0; i<5; i++) Serial.printf("%d: %d\n", i, coeffs[i]);
}

void loop() {
  delay(1000);
}

On the library as it stands, you get:

Starting
0: 131072
1: 0
2: 0
3: 0
4: 0

Which looks like the default: case values of a 'benign passthrough' function. With the simple patch of:

diff --git a/control_sgtl5000.cpp b/control_sgtl5000.cpp
index fe2c374..3f52b86 100644
--- a/control_sgtl5000.cpp
+++ b/control_sgtl5000.cpp
@@ -1000,6 +1000,7 @@ void calcBiquad(uint8_t filtertype, float fC, float dB_Gain, float Q, uint32_t q
     a0 = (A+1.0F) - ((A-1.0F)*cosw) + (beta*sinw);
     a1 = -2.0F * ((A-1.0F) - ((A+1.0F)*cosw));
     a2 = -((A+1.0F) - ((A-1.0F)*cosw) - (beta*sinw));
+  break;
   default:
     b0 = 0.5;
     b1 = 0.0;

you get:

Starting
0: 262144
1: -471615
2: 214299
3: 471616
4: -214298

which is looking much more like what you'd probably expect.

I'll send the trivial PR in a moment.

[robsoles]

Well, that's embarrassing - I have my old repository on github reflecting the omitted break; I would swear I saw 'rational' results for testing each filter back then but either I fixed it for my ex-employers and forgot to even fix my repository, let alone make a pull request if my memory of seeing them all work is not false. Cheers for clearing that up, did you include coeff flips for those two coefficients? It may as well be done in the final assignment of each in the calc_biquad routine itself imho - when we switched to floating point processor we also switched codecs to have a better matching third input and I was very annoyed I made that codec work seamlessly with the audio library and never got to submit it to the library; it was a sh1t hot chip and I tidied and tightened everything up around 96K sampling of four channels for crisper audio and would have been stoked to share it :/ Anyway, does that "unsigned short AudioControlSGTL5000::autoVolumeControl(uint8_t maxGain, uint8_t lbiResponse, uint8_t hardLimit, float threshold, float attack, float decay)" routine of mine do anything rational? I was just casting my eye over my old repository copy [forgive me, I should look at the live copy to see if anyone has had to fix it or anything] but it looks at least a bit like I am inviting the compiler to do floating point calculations in an integer container [lol, my current boss is as much as paranoid about being really explicit about what you expect to happen] and in that return call to AudioControlSGTL5000::modify I wonder if I am inviting the compiler to shift 12 bits left in an 8 bit container in spite of that the target container specified in the 'modify' function model is a minimum of 16 bits even in your poorest systems - if I was writing that stuff in front of my current boss nothing would be shifted more bits than it itself is [regardless of target's format] and those floating point calc into 8 bit containers would be as explicit as possible; uint8_t target=(uint8_t)((float)source/this*that+whatever); gets the desired result in all compilers that claim to do C. I do not recall testing autoVolumeControl - it was never a candidate for use in our project... Rob.

…

On Tue, Dec 29, 2020 at 10:27 PM Graham Whaley ***@***.***> wrote: Hi @robsoles <https://github.com/robsoles> - thanks for responding. The bug isn't a problem with the calculations themselves, just with a missing break; at the end of the hishelf switch case, which means it falls through and wipes out all the calculated coeffs with the default (benign passthrough) ones... I knocked up a simple example ino to show it (as generally requested in the Issue template): #include <Audio.h> void setup() { int coeffs[5]; delay(1000); Serial.begin(115200); Serial.println("Starting"); calcBiquad(FILTER_HISHELF, 1000, 0.0, 0.707, 524288, AUDIO_SAMPLE_RATE_EXACT, coeffs); for( int i=0; i<5; i++) Serial.printf("%d: %d\n", i, coeffs[i]); } void loop() { delay(1000); } On the library as it stands, you get: Starting 0: 131072 1: 0 2: 0 3: 0 4: 0 Which looks like the default: case values of a 'benign passthrough' function. With the simple patch of: diff --git a/control_sgtl5000.cpp b/control_sgtl5000.cpp index fe2c374..3f52b86 100644--- a/control_sgtl5000.cpp+++ b/control_sgtl5000.cpp@@ -1000,6 +1000,7 @@ void calcBiquad(uint8_t filtertype, float fC, float dB_Gain, float Q, uint32_t q a0 = (A+1.0F) - ((A-1.0F)*cosw) + (beta*sinw); a1 = -2.0F * ((A-1.0F) - ((A+1.0F)*cosw)); a2 = -((A+1.0F) - ((A-1.0F)*cosw) - (beta*sinw));+ break; default: b0 = 0.5; b1 = 0.0; you get: Starting 0: 262144 1: -471615 2: 214299 3: 471616 4: -214298 which is looking much more like what you'd probably expect. I'll send the trivial PR in a moment. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#343 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABSHQUTSF5S2QSTKJ72PF3DSXG4JDANCNFSM4NCTEOOA> .

[grahamwhaley]

Well, that's embarrassing

Heh, happens to us all!

Cheers for clearing that up

np

did you include coeff flips for those two coefficients? It may as well be done in the final assignment of each in the calc_biquad routine itself imho

I didn't. I normally keep Issues and PRs small and related to one item (bug, issue) at a time, so if I pushed a PR it would be just for the -1 coefficient fix, and tied to Issue #372 . I was not intending to do a PR for that, given I don't really understand the math or exactly how biquad filters work. Sure, I could do it rote, but I've no simple way of checking I got it right to hand (but, one could probably be constructed with some tone generator components and some tone detect or fft blocks maybe).

We probably want some input from @PaulStoffregen on the best place to make the fix. Personally, factoring out the coeff calculations to a single piece of code might be nice, rather than having what look to be near duplicates for all the calculations in both the sg5k code and the biquad_filter code :-)

Anyway, does that "unsigned short AudioControlSGTL5000::autoVolumeControl(uint8_t maxGain, uint8_t lbiResponse, uint8_t hardLimit, float threshold, float attack, float decay)" routine of mine do anything rational? I was just casting my eye over my old repository copy [forgive me, I should look at the live copy to see if anyone has had to fix it or anything] but it looks at least a bit like I am inviting the compiler to do floating point calculations in an integer container

Heh, no idea!. Well, I do call it to set the AGC up in some other code (not yet up on github, but one day). I took the argument values from the one and only place I could find relevant numbers, from the example code. I did try to do a bit of research to find now to map what one might consider a 'normal' AGC for speech (low bandwidth speech - 3.5KHz amateur radio for instance) to the SGTL hardware settings, but failed to find much. So far those settings seem to do something and are not hurting my ears...

There is the open Issue #164 that looks vaguely related to that area btw.

[robsoles]

Good stuff Graham. I popped a closer scrutinise over autoVolumeControl and it is likely the compiler is nailing the intent of the code tho I may make it more explicit myself if I ever PR that repository again. With regards to issue #164, if I commented I would ask them why they think system maximum volume is a good target [or hard limit] where this is the maximum you can hit not just before you clip, in this case, but just as you clip for crying out loud - I strongly recommend negative figures for both of these, I would not send that routine anything closer to zero than -1 for target and -0.5 for hard limit; I nearly told Pensive on the forums but he sorted himself out; I may have been half dead when #!64 was raised. I'd say more rational figures would be -3 and -2. Using the signal as reference 'dB' is a scalar where positive figures indicate apply gain and negative figures apply attenuation, when referring to a system's 'power' 'dB' is its scale where 0dB is the [hard] limit of said power's scale but in effect it remains scalar to any input signal but only in attenuation. So, fail to clip your input buffer all the way to the DAC and 0dB will not clip; keep your input buffer below -#dB and you can apply positive #dB gain figure to volume and still not clip. Sorry if this has all been relatively elementary to you. As to calc_biquad, funnily enough it came to me while I was on the verge of sleep and [this could be wrong but] I think the SGTL5000's inbuilt biquad stuff prefers the flipped signs on those two coefficients and that is what made me leave it in the first place - as the quantization figure is clearly different for that scenario sign flipping could be handled in calc_biquad based on what value is passed for that. I almost started understanding how the blasted bi-quadratic equations were working while I wrote the ASM to make what I wanted it to do to be as efficiently as possible, I could almost write a variant of the biquad processing block off the top of my head in C/++ today; but I would not as I always err in favor of refreshers for greater certainty... For a laugh, I provide that I am currently working on a project wherein I intend to bore an ~11mm hole into the side of a largish transformer [240|+15_0_-15/24/12 or something like that] and set an old teensy 3.1 up to control the motor turning the shaft with the N52 magnet on the end inserted with correct orientation into the hole in the transformer with intention to have my code slowly rise from 1Hz to 1000Hz [revolutions of magnet per second, I have to tag the shaft and count visually or I will add drag/imbalance etc] whilst reading the largest winding via a suitable voltage divider while that winding is connected to an eight ohm coil I wound on purpose just to have a relatively harmless 8 ohm load I could abuse the sh1t out of without caring too much - yep, I am going hunting for the resonant [hopefully 'take off'] frequency; once I've mapped 1Hz to 1KHz if I find a sweet spot I will try to set up to investigate harmonics of that sweet spot to see if I can get nuttiness to occur - of course, at some point I wish to fully bridge rectify the 24V winding, tank it in a nice big cap [48V 6.8mF standing by, lol] and then feed it to the drive circuit's input voltage point via a diode with something ridiculous coming off that 240V winding that I can use to tease Physicists and Electronics experts alike - even if I only just get it to self-sustain by using the 240V winding instead of the 24V one I will still have plenty to show them; I may have to slit the side of the transformer I bore the hole in for best effect but I wish to see it intact first as I suspect that the field will spin nicely without the split and may suffer less drag just because there are no gaps around the magnet block. Alternatively, it will be 'the right sized gap causes the greater slap yielding higher output' but I would prefer two identical transformers cut to each idea before attempting to jump to conclusions about this - I was planning a 3 phase construction of my own when I realised the good old fashioned transformer was what I was going to emulate for each of my phases anyway - I have already established I can get interesting voltages with a rod of iron bent around to point back at itself and wound ~120 turns ~2ohms. N36 magnet put 2V [bridge rectified and loaded ~100mA into an LED array] out of that at merely 300Hz - of course I was pouring about 4V @400mA+ to achieve this but my iron rod was clearly giving the magnet a 'slap' every half turn as it hardly presented a 'balanced' area for the magnet to turn in. Well, I laughed, perhaps you are just facepalming? :D Best, Rob.

…

On Wed, Dec 30, 2020 at 3:57 AM Graham Whaley ***@***.***> wrote: Well, that's embarrassing Heh, happens to us all! Cheers for clearing that up np did you include coeff flips for those two coefficients? It may as well be done in the final assignment of each in the calc_biquad routine itself imho I didn't. I normally keep Issues and PRs small and related to one item (bug, issue) at a time, so if I pushed a PR it would be just for the -1 coefficient fix, and tied to Issue #372 <#372> . I was not intending to do a PR for that, given I don't really understand the math or exactly how biquad filters work. Sure, I could do it rote, but I've no simple way of checking I got it right to hand (but, one could probably be constructed with some tone generator components and some tone detect or fft blocks maybe). We probably want some input from @PaulStoffregen <https://github.com/PaulStoffregen> on the best place to make the fix. Personally, factoring out the coeff calculations to a single piece of code might be nice, rather than having what look to be near duplicates for all the calculations in both the sg5k code and the biquad_filter code :-) Anyway, does that "unsigned short AudioControlSGTL5000::autoVolumeControl(uint8_t maxGain, uint8_t lbiResponse, uint8_t hardLimit, float threshold, float attack, float decay)" routine of mine do anything rational? I was just casting my eye over my old repository copy [forgive me, I should look at the live copy to see if anyone has had to fix it or anything] but it looks at least a bit like I am inviting the compiler to do floating point calculations in an integer container Heh, no idea!. Well, I do call it to set the AGC up in some other code (not yet up on github, but one day). I took the argument values from the one and only place I could find relevant numbers, from the example code <https://github.com/PaulStoffregen/Audio/blob/master/examples/HardwareTesting/SGTL5000/dap_avc_agc/dap_avc_agc.ino#L42>. I did try to do a bit of research to find now to map what one might consider a 'normal' AGC for speech (low bandwidth speech - 3.5KHz amateur radio for instance) to the SGTL hardware settings, but failed to find much. So far those settings seem to do something and are not hurting my ears... There is the open Issue #164 <#164> that looks vaguely related to that area btw. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#343 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABSHQUQJ67DBBA2A55IJARLSXIDAJANCNFSM4NCTEOOA> .