LowpassFourPole UnitTest has inaccurate Impulse Response

[tap]

Might be due to some roundoff in Processing?

@nwolek do you have the code you used to generate the reference response?

[nwolek]

Looks like Processing has spotty support for doubles: "A double can be used similarly to a float, but can have a greater magnitude because it's a 64-bit number. Processing functions don't use this datatype, so while they work in the language, you'll usually have to convert to a float using the (float) syntax before passing into a function."

That's a shame, because it was so easy to write the following code in Processing. I will look for another solution, because I think we need to be in the habit of building our targets outside of jamoma2.

double kPi  = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068;
double kTwoPi  = kPi * 2.0;
double in1 = 0;
double in2 = 0;
double in3 = 0;
double in4 = 0;
double out1 = 0;
double out2 = 0;
double out3 = 0;
double out4 = 0;
double frequency = 1000.0;
double sampleRate = 44100.0;
double tenthResonance = 1.0;
// linFreq 0-1 maps to 0-half the sampling rate
//double radians = (5000. * kTwoPi) / 44100.;
double fNormalizedToHalfNyquist = frequency * 4 / sampleRate;
void setup() {
  println("fNormalizedToHalfNyquist is " + fNormalizedToHalfNyquist);

  for (int i = 0; i < 64; i++)
  {
     double x, y;

      if (i==0) 
     {
        x = 1.0;
     } else {
        x = 0.0;
     }
     y = process(x, fNormalizedToHalfNyquist, tenthResonance);

     //println("x = " + x + " - y = " + y);
     println(y);
  }
}
double process(double input, double fc, double res)
{
  // source: http://musicdsp.org/archive.php?classid=3#26
  double f = fc * 1.4716;
  double fb = res * (1.0 - 0.4 * f * f);
  input -= out4 * fb;
  input *= 0.35013 * (f*f)*(f*f);
  out1 = input + 0.3 * in1 + (1 - f) * out1; // Pole 1
  in1  = input;
  out2 = out1 + 0.3 * in2 + (1 - f) * out2;  // Pole 2
  in2  = out1;
  out3 = out2 + 0.3 * in3 + (1 - f) * out3;  // Pole 3
  in3  = out2;
  out4 = out3 + 0.3 * in4 + (1 - f) * out4;  // Pole 4
  in4  = out3;
  return out4;
}

[lossius]

My suggestion would be that we do this in Matlab or Octave whenever we can.

[tap]

Agree with both of you. Octave is better than Matlab when possible because it is free, though if we need to use the Filter Toolbox or something then obviously we will need Matlab for those cases.

[nwolek]

I am not against Octave, I am just more fluent in Processing. Seems like the doubles issue will push us toward Octave though.

[nwolek]

Working today on changing the syntax to work in Octave. I have the following which writes output to the console. Next step will be to write output to a file.

clear
kPi  = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068;
kTwoPi  = kPi * 2.0;
in1 = 0;
in2 = 0;
in3 = 0;
in4 = 0;
out1 = 0;
out2 = 0;
out3 = 0;
out4 = 0;
frequency = 1000.0;
sampleRate = 44100.0;
tenthResonance = 1.0;
% linFreq 0-1 maps to 0-half the sampling rate
% radians = (5000. * kTwoPi) / 44100.;
fNormalizedToHalfNyquist = frequency * 4 / sampleRate;
for i = 0:64
    if (i==0)
        x = 1.0;
    else
        x = 0.0;
    endif
    %% source: http://musicdsp.org/archive.php?classid=3#26
    f = fNormalizedToHalfNyquist * 1.4716;
    fb = tenthResonance * (1.0 - 0.4 * f * f);
    in0 = x - (out4 * fb);
    in0 = in0 * (0.35013 * (f*f)*(f*f));
    out1 = in0 + 0.3 * in1 + (1 - f) * out1; % Pole 1
    in1  = in0;
    out2 = out1 + 0.3 * in2 + (1 - f) * out2;  % Pole 2
    in2  = out1;
    out3 = out2 + 0.3 * in3 + (1 - f) * out3;  % Pole 3
    in3  = out2;
    out4 = out3 + 0.3 * in4 + (1 - f) * out4;  % Pole 4
    in4  = out3;
    out4

endfor