1
$\begingroup$

I think the title says it all. Anyone could provide a solid answer to why it happens? Here is a sample code to replicate the issue.

    >> sos
sos =

    0.9293         0         0    1.0000         0         0
    1.0029   -1.9972    0.9944    1.0029   -1.9972    0.9944
    1.0139   -2.0190    1.0052    1.0138   -2.0190    1.0052
    0.9734   -1.9307    0.9576    0.9739   -1.9307    0.9570
    1.0274   -2.0193    0.9929    1.0271   -2.0193    0.9932
    1.0322   -1.9926    0.9647    1.0323   -1.9926    0.9645
    1.0721   -1.9903    0.9354    1.0716   -1.9903    0.9359
    1.1367   -1.9326    0.8640    1.1366   -1.9326    0.8641
    1.2812   -1.7360    0.7233    1.2800   -1.7360    0.7246
    1.5390   -0.9713    0.4036    1.5731   -0.9713    0.3696
    2.8934   -0.2280    0.4996    2.9959   -0.3522    0.5213

>> sos_norm

sos_norm =

    0.9293         0         0    1.0000         0         0
    1.0000   -1.9915    0.9915    1.0000   -1.9915    0.9915
    1.0001   -1.9915    0.9915    1.0000   -1.9915    0.9916
    0.9994   -1.9824    0.9832    1.0000   -1.9824    0.9826
    1.0003   -1.9660    0.9667    1.0000   -1.9660    0.9670
    0.9999   -1.9302    0.9345    1.0000   -1.9302    0.9344
    1.0005   -1.8574    0.8729    1.0000   -1.8574    0.8734
    1.0001   -1.7003    0.7601    1.0000   -1.7003    0.7602
    1.0010   -1.3563    0.5651    1.0000   -1.3563    0.5661
    0.9784   -0.6175    0.2566    1.0000   -0.6175    0.2349
    0.9658   -0.0761    0.1668    1.0000   -0.1176    0.1740

>> y1 = dfilt.df2sos(sos).filter(x);
>> y2 = dfilt.df2sos(sos_norm).filter(x);
>> sum(abs(y1 - y2))
ans =

   1.8913e-09

Upon the request I'm providing different code block to show the difference.

load("sos.mat");

sos = double(sos);
x = double(rand(1000,1));

sos_norm = zeros(size(sos));
for i = 1:size(sos, 1)
    sos_norm(i,:) = sos(i,:) ./ sos(i,4);
end

y1 = sosfilt(sos,x);
y2 = sosfilt(sos_norm,x);

disp("diff for sosfilt");disp(max( abs( y1 - y2 ) ) / max( abs( y1 ) ));

y3 = dfilt.df2sos(sos).filter(x);
y4 = dfilt.df2sos(sos_norm).filter(x);


disp("diff for dfilt filter");disp(max( abs( y3 - y4 ) ) / max( abs( y3 ) ));

Result:

diff for sosfilt
     0

diff for dfilt filter
   9.5757e-12
Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

The result sum(abs(y1 - y2)) doesn't tell you much, because its value depends on the length of the input signal. You could look at the value of

max( abs( y1 - y2 ) ) / max( abs( y1 ) )

and I'm pretty sure that you'll find out that this value is in the order of the round-off errors that you can expect using double precision computations. So what you see is just a consequence of doing computations with finite precision.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ I updated the code upon request. Please, let me know what do you think about the problem. $\endgroup$
    – Tərlan Əhəd
    Sep 18 at 11:41
  • 1
    $\begingroup$ The problem is it happens with dfilt but not with sosfilt. why? If problem is with precision then it should happen with sosfilt as well. $\endgroup$
    – Tərlan Əhəd
    Sep 18 at 11:44
  • 1
    $\begingroup$ @TərlanƏhəd: I don't know the exact internal computations in sosfilt vs. dfilt.df2sos. It could be that sosfilt internally normalizes the second-order coefficients, which would give you the result you see. But anyway, a relative error of 1e-12 is nothing to worry about, it is indeed a round-off error. $\endgroup$
    – Matt L.
    Sep 18 at 12:18
  • 1
    $\begingroup$ @TərlanƏhəd: I just confirmed that sosfilt does indeed normalize by $a_0$, that's why it doesn't make a difference whether or not you provide a normalized sos matrix or not. $\endgroup$
    – Matt L.
    Sep 18 at 12:21
1
$\begingroup$

First of all the difference is extremely small, it represents a numerical noise level of about -230dB, so it's indeed numerical noise. Keep in mind that you are using floating numbers, so all values are always "approximate" and never "precise".

The issue here is that an unnormalized second order section cannot be implemented directly: it has to be normalized before you can turn it into code. There are different methods of doing this and apparently sosfilt() and dfilt() use different ones. Obviously sosfilt() uses the same method that you do by normalizing the denominators only. Another popular choice is to normalize the numerators as well, i.e. normalize each section so that $b_0 = 1$ and put the cumulative gain into the first section. A normalization algorithm can also decide to re-order the section based on pole locations or adjust the grouping of poles and zeros which can have significant impact on the numerical noise.

Since neither sosfilt() or dfilt() have documented the exact method of normalization, we don't know for sure. Please keep in mind that none of these methods is inherently "right" or "wrong" they are just different and since there will be always numerical noise, the numerical noise will be different too.

Improve this answer
$\endgroup$

Your Answer

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

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