Wrong results when cascading biquads in sosfilter_c

[happyTonakai]

sosfilter_c gives wrong results when cascading multiple biquads. Some simple code to show:

import numpy as np
from pyfilterbank.sosfiltering import sosfilter_c
from scipy.signal import sosfilt

sos0 = np.array([
    [0.9994671940803528, -1.980954885482788, 0.9816480278968811, 1.0, -1.9809452295303345, 0.981124997138977],
    [0.9940603375434875, -1.8794806003570557, 0.89670330286026, 1.0, -1.8794806003570557, 0.8907635807991028]
])
sos1 = sos0.flatten()

x = np.random.random((16,))
z0 = np.zeros((sos0.shape[0], 2))
y0, z0 = sosfilt(sos0, x, zi=z0)

y1, z1 = sosfilter_c(x, sos1)

y2, _ = sosfilter_c(x, sos0[0])
y2, _ = sosfilter_c(y2, sos0[1])

y0 and y2 are identical except the numerical precision but y1 is incorrect.

[SiggiGue]

You should recheck your code (there is a wrong usage of the function) and read the documentation.

ZR Han @.***> schrieb am Fr., 6. Jan. 2023, 10:16:

sosfilter_c gives wrong results when cascading multiple biquads. Some simple code to show:

import numpy as npfrom pyfilterbank.sosfiltering import sosfilter_cfrom scipy.signal import sosfilt sos0 = np.array([ [0.9994671940803528, -1.980954885482788, 0.9816480278968811, 1.0, -1.9809452295303345, 0.981124997138977], [0.9940603375434875, -1.8794806003570557, 0.89670330286026, 1.0, -1.8794806003570557, 0.8907635807991028] ])sos1 = sos0.flatten() x = np.random.random((16,))z0 = np.zeros((sos0.shape[0], 2))y0, z0 = sosfilt(sos0, x, zi=z0) y1, z1 = sosfilterc(x, sos1) y2, = sosfilterc(x, sos0[0])y2, = sosfilter_c(y2, sos0[1])

y0 and y1 are identical except the numerical precision but y1 is incorrect.

— Reply to this email directly, view it on GitHub https://github.com/SiggiGue/pyfilterbank/issues/24, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAXQ6H7VV6QZFFIXDU5UVV3WQ7PIJANCNFSM6AAAAAATS37RWA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[happyTonakai]

@SiggiGue Thank you for response. Could you please specify where the wrong usage is?

I've read the documentation and the syntax of the function is

signal_out, states = sosfilter_c(signal_in, sos, states=None)

signal_in has a shape of (N, 0) and sos has a shape of (K*6, 0). states can be None.

In my code x has a shape of (16, ), sos1 has a shape of (2*6, ), and initial states is None. These usage are all met in my example code and I don't find anything goes wrong. I also tried to convert sos and x into np.float32 which doesn't work as well.

By the way, I think there's something wrong in your unit test code:

https://github.com/SiggiGue/pyfilterbank/blob/5d74af8d936d79dcb61f7f35a3120d79ea93480b/pyfilterbank/tests/test_sosfiltering.py#L10-L18

In the test code sos has a shape of (2, 6) and results in ksos=0 because len(sos)=2.

[happyTonakai]

You can check this out by replacing sf.sosfilter_double_c with sf.sosfilter_c in the ut https://github.com/SiggiGue/pyfilterbank/blob/87b9bb52cb6ecdbcbc0385f1c10ed040b22b49d1/pyfilterbank/tests/test_sosfiltering.py#L23 and you will see the unit test failed:

Traceback (most recent call last):
  File "~/ut.py", line 41, in test_implementation
    self.assertAlmostEqual(sop, soc, places=6)
AssertionError: 1119.841385113597 != 1121.028 within 6 places (1.1865689879655292 difference)

[SiggiGue]

If you cascade sos, you need too feed the intermediate states into the function as well to get the expected results for your comparison.

ZR Han @.***> schrieb am Do., 12. Jan. 2023, 03:40:

You can check this out by replacing sf.sosfilter_double_c with sf.sosfilter_c in the ut

https://github.com/SiggiGue/pyfilterbank/blob/87b9bb52cb6ecdbcbc0385f1c10ed040b22b49d1/pyfilterbank/tests/test_sosfiltering.py#L23 and you will see the unit test failed:

Traceback (most recent call last): File "~/ut.py", line 41, in test_implementation self.assertAlmostEqual(sop, soc, places=6) AssertionError: 1119.841385113597 != 1121.028 within 6 places (1.1865689879655292 difference)

— Reply to this email directly, view it on GitHub https://github.com/SiggiGue/pyfilterbank/issues/24#issuecomment-1379739305, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAXQ6H4657QOED2Z6MT4BWTWR5VIRANCNFSM6AAAAAATS37RWA . You are receiving this because you were mentioned.Message ID: @.***>

[happyTonakai]

If you cascade sos, you need too feed the intermediate states into the function as well to get the expected results for your comparison.

This doesn't make sense to me. The documentation states that initial state can be None, and even I change this line

y1, z1 = sosfilter_c(x, sos1)

to

y1, z1 = sosfilter_c(x, sos1, states=np.zeros((2 * 2, )))

the result is still incorrect. If you meant that

y2, _ = sosfilter_c(x, sos0[0])
y2, _ = sosfilter_c(y2, sos0[1])

should be

y2, state = sosfilter_c(x, sos0[0])
y2, _ = sosfilter_c(y2, sos0[1], state)

I don't think this is true. Cascaded biquads don't share internal states. y2 gives the same result as scipy.signal.sosfilt.

[happyTonakai]

Have a look at this code:

import numpy as np
from pyfilterbank.sosfiltering import sosfilter_c, sosfilter_double_c
from scipy.signal import sosfilt

sos = np.array([
    [0.9994671940803528, -1.980954885482788, 0.9816480278968811, 1.0, -1.9809452295303345, 0.981124997138977],
    [0.9940603375434875, -1.8794806003570557, 0.89670330286026, 1.0, -1.8794806003570557, 0.8907635807991028]
], dtype=np.float32)

x = np.random.random((16,)).astype(np.float32)

y0, z0 = sosfilt(sos, x, zi=np.zeros((sos.shape[0], 2)))
y1, z1 = sosfilter_c(x, sos.flatten(), states=np.zeros((sos.shape[0] * 2, )))
y2, z2 = sosfilter_double_c(x, sos.flatten(), states=np.zeros((sos.shape[0] * 2, )))

print(np.mean((y0 - y1)**2))
print(np.mean((y0 - y2)**2))

The only difference between y1 and y2 is sosfilter_c and sosfilter_double_c, and the output is quite different.

0.0013462692831356475
3.282670120428347e-29

[SiggiGue]

Sorry, you are totally right with the states, I was a bit short in time, answering to your issue.

Maybe you can find out, what the problem is? Because I have no spare time for this repo at the moment.

The plan is to remove the c part totally since scipy has a sosfilt function by this time.

Best regards and thank you for reporting.

ZR Han @.***> schrieb am Do., 12. Jan. 2023, 08:56:

Have a look at this code:

import numpy as npfrom pyfilterbank.sosfiltering import sosfilter_c, sosfilter_double_cfrom scipy.signal import sosfilt sos = np.array([ [0.9994671940803528, -1.980954885482788, 0.9816480278968811, 1.0, -1.9809452295303345, 0.981124997138977], [0.9940603375434875, -1.8794806003570557, 0.89670330286026, 1.0, -1.8794806003570557, 0.8907635807991028] ], dtype=np.float32) x = np.random.random((16,)).astype(np.float32) y0, z0 = sosfilt(sos, x, zi=np.zeros((sos.shape[0], 2)))y1, z1 = sosfilter_c(x, sos.flatten(), states=np.zeros((sos.shape[0] 2, )))y2, z2 = sosfilter_double_c(x, sos.flatten(), states=np.zeros((sos.shape[0] 2, ))) print(np.mean((y0 - y1)2))print(np.mean((y0 - y2)2))

The only difference between y1 and y2 is sosfilter_c and sosfilter_double_c, and the output is quite different.

0.0013462692831356475 3.282670120428347e-29

— Reply to this email directly, view it on GitHub https://github.com/SiggiGue/pyfilterbank/issues/24#issuecomment-1379938263, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAXQ6HYSD2ZKMQQFXUKNSG3WR62LXANCNFSM6AAAAAATS37RWA . You are receiving this because you were mentioned.Message ID: @.***>

[SiggiGue]

Have you tried what happens if you feed the scipy sosfilt function with float32 data cascaded inside and outside? Maybe this is just a numerical issue.

Siegfried Gündert @.***> schrieb am Do., 12. Jan. 2023, 09:45:

Sorry, you are totally right with the states, I was a bit short in time, answering to your issue.

Maybe you can find out, what the problem is? Because I have no spare time for this repo at the moment.

The plan is to remove the c part totally since scipy has a sosfilt function by this time.

Best regards and thank you for reporting.

ZR Han @.***> schrieb am Do., 12. Jan. 2023, 08:56:

Have a look at this code:

import numpy as npfrom pyfilterbank.sosfiltering import sosfilter_c, sosfilter_double_cfrom scipy.signal import sosfilt sos = np.array([ [0.9994671940803528, -1.980954885482788, 0.9816480278968811, 1.0, -1.9809452295303345, 0.981124997138977], [0.9940603375434875, -1.8794806003570557, 0.89670330286026, 1.0, -1.8794806003570557, 0.8907635807991028] ], dtype=np.float32) x = np.random.random((16,)).astype(np.float32) y0, z0 = sosfilt(sos, x, zi=np.zeros((sos.shape[0], 2)))y1, z1 = sosfilter_c(x, sos.flatten(), states=np.zeros((sos.shape[0] 2, )))y2, z2 = sosfilter_double_c(x, sos.flatten(), states=np.zeros((sos.shape[0] 2, ))) print(np.mean((y0 - y1)2))print(np.mean((y0 - y2)2))

The only difference between y1 and y2 is sosfilter_c and sosfilter_double_c, and the output is quite different.

0.0013462692831356475 3.282670120428347e-29

— Reply to this email directly, view it on GitHub https://github.com/SiggiGue/pyfilterbank/issues/24#issuecomment-1379938263, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAXQ6HYSD2ZKMQQFXUKNSG3WR62LXANCNFSM6AAAAAATS37RWA . You are receiving this because you were mentioned.Message ID: @.***>

[happyTonakai]

As you can see the output is so different that we can't say it's related to the numeric issue. The states of the first stage of sos is correct but the second one is not. I've checked the source code and all I found is

https://github.com/SiggiGue/pyfilterbank/blob/87b9bb52cb6ecdbcbc0385f1c10ed040b22b49d1/pyfilterbank/sosfiltering.py#L109-L110

can be changed to

 if isinstance(states, type(None)): 
     states = np.zeros(ksos*2).astype(np.float32) 

but this doesn't help because you convert it to np.float32 in the following lines

https://github.com/SiggiGue/pyfilterbank/blob/87b9bb52cb6ecdbcbc0385f1c10ed040b22b49d1/pyfilterbank/sosfiltering.py#L112-L113

C code seems correct too. I can't tell what's going wrong.

[SiggiGue]

Thank you for having a look into that. Line 111 is a bit misleading withe the double cast and I should have changed it to float32 but as you said, this won't make a difference.

I would be very interested in the reson of this issue. Anyway, applying the YAGNI principle will remove the sosfilt.c stuff from this repo, as mentioned above.

ZR Han @.***> schrieb am Do., 12. Jan. 2023, 10:30:

[image: image] https://user-images.githubusercontent.com/21305646/212028051-2830982d-3e86-40b0-b4cf-fb5efa33b988.png

As you can see the output is so different that we can't say it's related to the numeric issue. The states of the first stage of sos is correct but the second one is not. I've checked the source code and all I found is

https://github.com/SiggiGue/pyfilterbank/blob/87b9bb52cb6ecdbcbc0385f1c10ed040b22b49d1/pyfilterbank/sosfiltering.py#L109-L110

can be changed to

if isinstance(states, type(None)): states = np.zeros(ksos*2).astype(np.float32)

but this doesn't help because you convert it to np.float32 in the following lines

https://github.com/SiggiGue/pyfilterbank/blob/87b9bb52cb6ecdbcbc0385f1c10ed040b22b49d1/pyfilterbank/sosfiltering.py#L112-L113

C code seems correct too. I can't tell what's going wrong.

— Reply to this email directly, view it on GitHub https://github.com/SiggiGue/pyfilterbank/issues/24#issuecomment-1380040802, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAXQ6H6H5CR6HNIKY7GDIHLWR7FKDANCNFSM6AAAAAATS37RWA . You are receiving this because you modified the open/close state.Message ID: @.***>