0
$\begingroup$

I'm trying to create a LPF, using language JAVA. My first aim is to understand how filters design is working and the second one is to create a "home made" equalizer with "home made" filters.
To help me, I'm using EarLevel website for the formulas and to check my calculations results for the feedforward / feedback coefficients ; and the introduction to digital filters by Julius O. Smith.
My input have is sampled at 44100 Hz and my Q=0.7071 (an approximation of 1/sqrt(2)) and I'm using the general formula :

$$y(n) = f_0 x(n) + f_1 x(n-1) + f_2 x(n-2) - b_1 y(n-1) - b_2 y(n-2)$$

with $f_0$ being the zeroth feed-forward coefficient, and so on, and $b_1$ being the first feed-back coefficient.
Edit : The coefficients are stored in a type "double".

In a first time, I set the center frequency at 10 000 Hz and it seems working. Then, I tried to lower the center frequency. When I reached around 7000 Hz, my output signal started to diverge after less than 30 "computations".

Here are my feedforward / feedback coefficients calculated for a 44100 Hz sampling signal with a central frenquency at 2000 Hz and Q = 0.7071 :

\begin{align} f_0 &= 0.016819123361395717\\ f_1 &= 0.033638246722791434\\ f_2 &= 0.016819123361395717\\ b_1 &= -1.601089848140876\\ b_2 &= 0.668366341586459 \end{align}

There are the same values thand those given by the site calculator so, I suppose they are correct.

But here is the input / output result after 20 signals sent to the filter :
Edit : The input datas are integers, representing the 16 bits input audio signal (from -32000 to +32000). The audio format used is a PCM_signed signal, sampled at 44100 Hz, 16 bits, 2 channels. I made test by comparison to check if the transformation PCM coded => integer => PCM coded was ok and, yes, it is ; so, the problem don't seems to come from there.
The output datas are stored in a double to keep the precision after the unit digit. Of course, I round them to an integer before sending them to the sound card.

id = 20 - input = -2997 output= -4511.426234017688
id = 21 - input = -3022 output= 97.52810634470279
id = 22 - input = -3506 output= -6592.151428768328
id = 23 - input = -3176 output= 2949.256088740449
id = 24 - input = -3272 output= -10840.64632338293
id = 25 - input = -3287 output= 8909.22977931816
id = 26 - input = -2555 output= -19536.606226891832
id = 27 - input = -3123 output= 21316.24417688955
id = 28 - input = -1890 output= -37414.30445881503
id = 29 - input = -2821 output= 47023.28302742664
id = 30 - input = -1497 output= -74302.58304596202
id = 31 - input = -2557 output= 100154.15317482647
(...)
id = 45 - input = -3401 output= 1.6528078507359352E7
id = 46 - input = -2367 output= -2.3763232415918328E7
id = 47 - input = -3865 output= 3.415084117779622E7
(...)
id = 2047 - input = 1971 output= Infinity
id = 2048 - input = 4435 output= -Infinity
id = 2049 - input = 2648 output= Infinity

What's wrong ? Why does that signal diverge ? And why around 7000 Hz ? A LPF should work for much lower frequencies if I want to redirect it to a subwoofer for example, isn't it ?

Thank you for your help...

Edit2 : Here are the methods I'm using to find the coefficients and the method I'm using to calculate the output. 

protected int centerFrequency, samplingFrequency;
protected double settingFactor;
protected double feedforward0, feedforward1, feedforward2, feedback1, feedback2;
protected double normalized, coefK;
private int test=0;

public LowPassFilter(int samplingFrequency, int centerFrequency, double qualityFactor) {
        super(samplingFrequency, centerFrequency, qualityFactor);
        coefK=Math.tan(Math.PI*centerFrequency/samplingFrequency);
        computeGains(settingFactor);
    }

    @Override
    protected void computeGains(double qualityFactor) {
        normalized= 1 / (1 + coefK / qualityFactor + coefK * coefK);
        feedforward0 = coefK * coefK * normalized;
        feedforward1 = 2 * feedforward0;
        feedforward2 = feedforward0;
        feedback1 = 2 * (coefK * coefK - 1) * normalized;
        feedback2 = (1 - coefK / qualityFactor + coefK * coefK) * normalized;
        System.out.println("LPF.computeGains : "+feedforward0+" ; "+feedforward1+" ; "+feedforward2+" ; "+feedback1+" ; "+feedback2);
    }

And here is the routine : 

protected LinkedList<Integer> inputValues2 = new LinkedList<Integer>();
protected LinkedList<Double> outputValues2 = new LinkedList<Double>();

public double process(int inputValue) {
        // On renseigne la valeur d'entrée dans le tableau des valeurs d'entrées
        inputValues2.removeFirst();
        inputValues2.add(inputValue);
        double outputValue= feedforward0*inputValues2.get(2) + feedforward1*inputValues2.get(1) + feedforward2*inputValues2.get(0) - feedback1*outputValues2.get(1) - feedback2*outputValues2.get(0);
        outputValues2.removeFirst();
        outputValues2.add(outputValue);
        if ((test>0 && test<100) || (test>2030 && test<2050)) {
            System.out.println("BiquadFilter.process - id = "+test+" - input = "+inputValues2.get(0)+" output= "+outputValues2.get(0)         );
        }
        test++;
        return outputValue;
    }
Improve this question
$\endgroup$
14
  • $\begingroup$ It's really cool that you're using a biquad, because it's a very mighty architecture for efficient implementation and hence widely used in audio processing, but it would be helpful to know what you're planning to do with that filter afterwards, to understand what you need. $\endgroup$
    – Marcus Müller
    Commented Dec 13, 2018 at 20:02
  • $\begingroup$ also, I might be missing something, but what do your input and output values represent? $\endgroup$
    – Marcus Müller
    Commented Dec 13, 2018 at 20:02
  • $\begingroup$ Thank you for your reply Marcus. I edited my post to explain my "plans". My input values are the 16 bit audio signal (from -32000 to +32000) and my output value is the audio signal after the filtering. I kept is with digits after the unity (and not converted to integer) to have a better precision. $\endgroup$
    – Dr_Click
    Commented Dec 13, 2018 at 20:21
  • $\begingroup$ 1 - What language are you using? 2 - What size are your variables. Inputs, outputs, coefficients, states. It is imperative to know. $\endgroup$
    – Ben
    Commented Dec 13, 2018 at 20:26
  • $\begingroup$ but these are merely 47 samples then, @Dr_Click, how can you derive a low pass (or not low pass) behaviour from these? Where do your 7000 Hz appear here? I'm really a bit confused. $\endgroup$
    – Marcus Müller
    Commented Dec 13, 2018 at 20:38

1 Answer 1

Reset to default
3
$\begingroup$

Let's try this in Octave/Matlab:

a = [1,-1.601089848140876,0.668366341586459];
b = [0.016819123361395717,0.033638246722791434,0.016819123361395717];
r = roots(a);
abs(r)
ans =

   0.81754
   0.81754

Since the roots are all inside the unit circle, your filter is stable.

The magnitude response looks like it should for a Butterworth filter. However, as you've noticed yourself, for $f_s=44100$ Hz, the $3$-dB cut-off frequency is $2000$ Hz:

enter image description here

I conclude that the problem does not come from the filter coefficients. You should show us your filtering routine, because if there is no silly input/output problem, then the problem must be in the filtering routine. Even though you say you just implemented the recursion shown in your question, believe me that I've seen lots of buggy filter implementations. As long as you don't show your implementation, nobody can really help you find the problem.

It should also be easy enough to use any free software package (Octave/SciPy/etc.) with the same input and the same filter coefficients, and compare the outputs. This will convince you (if still necessary) that the problem is not caused by the coefficients.

You could perform a first simple check of your routine by computing the impulse response, using an impulse as the input signal (i.e., a 1 followed by zeros). In Octave/Matlab you can find the impulse response using the command impz (given the coefficients a and b from above):

h = impz(b,a);

The figure below shows the first $60$ values of the impulse response. Clearly, the filter is stable.

enter image description here

The first $15$ values of the impulse response are

   0.0168191
   0.0605672
   0.1025513
   0.1237127
   0.1295334
   0.1247092
   0.1130948
   0.0977236
   0.0808754
   0.0641737
   0.0486934
   0.0350710
   0.0236068
   0.0143563
   0.0072078
Improve this answer
$\endgroup$
8
  • $\begingroup$ Thank you for that first answer. I will try Octave / Matlab this week-end. I edited my question and added the routine I'm using to calculate my output. Do you see something wrong in it ? That's often true that I don't see my own errors while they would be obvious for a third-part reader... $\endgroup$
    – Dr_Click
    Commented Dec 15, 2018 at 9:57
  • $\begingroup$ @Dr_Click: You should use an impulse as input (a single value of 1 and then only zeros), and see what the output is. It should be the filter's impulse response, and it is easy to check if it's correct. $\endgroup$
    – Matt L.
    Commented Dec 15, 2018 at 17:22
  • $\begingroup$ I wrote a method to get the impulse response of my filter. But, how to know if the result is correct ? Is there a special site ? Or do I have to compare in MathLab / Octave ? $\endgroup$
    – Dr_Click
    Commented Dec 16, 2018 at 9:38
  • $\begingroup$ @Dr_Click: I added some info about the impulse response to my answer; use those to check your result. $\endgroup$
    – Matt L.
    Commented Dec 16, 2018 at 11:45
  • $\begingroup$ Ok, I don't get the same impulse result. You are right, there is (more than) probably a mistake in my routine. When you read my code, did you fond something weird or a mistake ? $\endgroup$
    – Dr_Click
    Commented Dec 16, 2018 at 22:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.