bugGNU Octave - Bugs: bug #63440, [octave forge] (signal) sosfilt -...

 
 

bug #63440: [octave forge] (signal) sosfilt - NaN

Submitter:  Renato <yamane>
Submitted:  Tue 29 Nov 2022 12:18:17 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  lostbard
Originator Name:  Open/Closed:  * Closed
Release:  * other Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 17 May 2023 05:52:11 PM UTC, comment #14: 

v1.4.4 released - closing ticket

John Donoghue <lostbard>
Group Member
Tue 25 Apr 2023 11:57:49 AM UTC, comment #13: 
John Donoghue <lostbard>
Group Member
Mon 27 Feb 2023 05:04:23 PM UTC, comment #12: 

@ Mike Miller

Did you see this bug report and the proposed solution (comment #3)?

I'm convinced that the modified functions zp2sos and sos2tf described in comment #3 (appended files 54164 to 54167) are an appreciable improvement and need to be included as soon as possible in the next update of octave and the signal processing package.

Thank you

Charles Praplan <charprap>
Wed 11 Jan 2023 05:14:51 PM UTC, comment #11: 

Hi Renato,
Unfortunately the function you want is not written. So you have to manage without it.
Independently of this function, just run the following code to plot the spectrum
before and after your while loop.
Although the while loop does its job of changing the crest factor, you will see that is
adding some LF (white) noise to your signal.



close all
clear
clc

pkg load ltfat

sampling_rate = 44100;
length_fact = 30;
hpf = 40;
lpf = 400;
crest_factor = 6;
filter_order = 4;
rounding = 4;

typenoise = noise((length_fact/2)*sampling_rate, 1, 'pink');
%load typenoise.mat

[z, p, k] = butter(filter_order, [hpf/(sampling_rate/2), lpf/(sampling_rate/2)]);
sos = zp2sos (z, p, k);
filtered = sosfilt(sos, typenoise);
normalized = filtered / (rms(filtered) / 10^(-crest_factor/20));
normalized_before_while_loop=normalized;

while crestfactor(normalized) > crest_factor
normalized(normalized > 0) = ((normalized(normalized > 0).^rounding) ./ ((normalized(normalized > 0).^rounding) + 1)).^(1/rounding);
normalized(normalized < 0) = (((normalized(normalized < 0).^rounding) ./ ((normalized(normalized < 0).^rounding) + 1)).^(1/rounding)) *-1;
rounding = rounding^2;
  if rounding > 100
    rounding = 100;
  end
  normalized = normalized / (rms(normalized) / 10^(-crest_factor/20));
end;

second_half = normalized * -1;
normalized = [[normalized]; [second_half]];

hist(normalized, 61);
xlabel('Amplitude');
ylabel('Samples');
grid on;

[spectra,freq] = pwelch(normalized_before_while_loop,hann(length(normalized_before_while_loop)),0,length(normalized_before_while_loop),sampling_rate);
figure
loglog(freq,spectra)
[spectra,freq] = pwelch(normalized,hann(length(normalized)),0,length(normalized),sampling_rate);
hold on
loglog(freq,spectra)
grid on
legend('spectrum before while loop', 'spectrum after while loop')
xlabel('[Hz]')
ylabel('Amplitude')


Charles Praplan <charprap>
Tue 10 Jan 2023 03:34:20 PM UTC, comment #10: 

comment #8:

> ...I don't know what your while loop does, but it is likely it is adding some white noise to your signal.
> I hope this helps.
>
> By the way, why don't you use octave to plot the spectrum?


Hi Charles,

The loop is changing the Crest Factor from 12dB to 6dB, that means the peak now is twice the average (peak = 1, average = 0.5), similar to what is explanined on Eq. 1. Page 3 of this paper (https://xlrtechs.com/dbkeele.com/PDF/Keele%20(2001-09%20AES%20%20Preprint)%20-%20Dev.%20Test%20Signals%20EIA-426-B.PDF)

Regarding the graph generated on Octave, it should be much better for me if the POCTAVE funcion was available (https://savannah.gnu.org/bugs/?63442)

Renato <yamane>
Sun 08 Jan 2023 03:02:52 PM UTC, comment #9: 

@ Mike Miller

Did you see this bug report and the proposed solution (comment #3)?

I'm convinced that the modified functions zp2sos and sos2tf described in comment #3 (appended files 54164 to 54167) are an appreciable improvement and need to be included as soon as possible in the next update of the signal processing package.

Thank you

Charles Praplan <charprap>
Wed 04 Jan 2023 05:48:47 PM UTC, comment #8: 

If you perform the same operation on the signal "normalized" BEFORE your while loop, and with the signal amplitude reduced to +-1 (because of fct audiowrite) by dividing "normalized" by 3 for example, you will see that the filter does its job!

I don't know what your while loop does, but it is likely it is adding some white noise to your signal.
I hope this helps.

By the way, why don't you use octave to plot the spectrum?

Charles Praplan <charprap>
Wed 04 Jan 2023 02:16:13 PM UTC, comment #7: 

Hi Charles,

Could you please confirm on your Matlab if this abnormal high-pass filter is also happening there?

At the end of the code on comment #4, please add the aditional line:


audiowrite ('signal.wav', normalized, sampling_rate);


And then you can open it on Audacity and go to menu Analyze --> Plot Spectrum

Looks like the low-pass filter is working ok, but the hi-pass filter is working just at a certain limit.


Renato <yamane>
Wed 04 Jan 2023 01:24:04 PM UTC, comment #6: 

Hi Charles,

I don't know exactly why, but today I was unable to reproduce the same strange histogram from yesterday, and unfortunatly I didn't save the output to look deeper.

So, I agree with you the proposed change is working okay.

Thanks for your effort to solve this bug.

Renato <yamane>
Wed 04 Jan 2023 09:58:38 AM UTC, comment #5: 

Hi Renato,

I made some comparisons between octave 7.3.0 (with latest zp2sos), octave 6.4.0 and matlab.

I checked the following variables and I could not see any differences (except some rounding at 10e-9, and ordering of values):

- z, p, k
- sos

If this is ok, the problem cannot be in zp2sos!

However, I also compared the variables "filtered" and "normalized" (before your while loop) using a signal "typenoise" that I saved in a file (else it is different at each execution and it is not available in Matlab).

Here again I could not see any problem.

So I'm pretty sure that your problem is not linked to the zp2sos function.

Can you make the same comparisons on your side, of course using the same data for all your tests (i.e. saving typenoise in a file)?




Charles Praplan <charprap>
Tue 03 Jan 2023 09:34:13 AM UTC, comment #4: 

Hi Charles,

Sorry for my late reply. I was on holiday.

I replaced the files zp2sos and sos2tf with your suggested ones and I noticed the histogram of the audiofile generated is very different between the original and the proposed changes.

Please test with this code:


sampling_rate = 44100;
length = 30;
hpf = 40;
lpf = 400;
crest_factor = 6;
filter_order = 4;
rounding = 4;

typenoise = noise((length/2)*sampling_rate, 1, 'pink');
[z, p, k] = butter(filter_order, [hpf/(sampling_rate/2), lpf/(sampling_rate/2)]);
sos = zp2sos (z, p, k);
filtered = sosfilt(sos, typenoise);
normalized = filtered / (rms(filtered) / 10^(-crest_factor/20));

while crestfactor(normalized) > crest_factor
normalized(normalized > 0) = ((normalized(normalized > 0).^rounding) ./ ((normalized(normalized > 0).^rounding) + 1)).^(1/rounding);
normalized(normalized < 0) = (((normalized(normalized < 0).^rounding) ./ ((normalized(normalized < 0).^rounding) + 1)).^(1/rounding)) *-1;
rounding = rounding^2;
  if rounding > 100
    rounding = 100;
  endif
  normalized = normalized / (rms(normalized) / 10^(-crest_factor/20));
endwhile;

second_half = normalized * -1;
normalized = [[normalized]; [second_half]];

hist(normalized, 61);
xlabel('Amplitude');
ylabel('Samples');
grid on;



Renato <yamane>
Wed 28 Dec 2022 10:57:36 AM UTC, comment #3: 

After some further testing, I propose following solution for the function

zp2sos:

  Problem:
  With current implementation of zp2sos, sos(i, 4) may be zero and this is
  unfortunately not supported by the fuction sosfilt.

  Solution:
   - The real poles and zeros must be combined (as the complex conjugates ones)
   - Leading zeros present in numerator and denominator are removed
     [0, 1, 1, 0, 1, 2] -> [1, 1, 0, 1, 2, 0]


Additionally, a bug has been found in

sos2tf:

  Problem:
  If several trailing zeros were present in numerator and denominator, only one
  of them were removed in numerator (index for A was not incremented)

  Correction:
  - same index for A as for B (line 87) because index for A was not incremented


Updated functions and corresponding diff files are appended!




(file #54164, file #54165, file #54166, file #54167)

Charles Praplan <charprap>
Sat 10 Dec 2022 10:36:11 AM UTC, comment #2: 

The case where 2 poles or zeros are real and equal, should probably be treated as a pair of complex conjugate.

Charles Praplan <charprap>
Sat 10 Dec 2022 09:44:15 AM UTC, comment #1: 

In https://savannah.gnu.org/bugs/?func=detailitem&item_id=60597, I proposed updates for several functions handling poles and zeroes to corrected several problems. This includes the function zp2sos.

During the development the following question arized : should the real poles and zeros be combined as the complex conjugates ones?

When rewriting the function zp2sos, the first idea was to have an input parameter, DoNotCombineReal, in the function \octave-7.3.0-w64\mingw64\share\octave\packages\signal-1.4.3\zp2sos.m to let the user decide.

Finally this parameter has been fixed to 1 inside the function (line 78).

Let examine the results in both cases

In your example,

Without combining the real poles and zeros we get (with your data):
  sos =
          0   0.0008   0.0008   1.0000  -1.9263   0.9296
          0   1.0000   1.0000   1.0000  -1.9914   0.9914
          0   1.0000  -1.0000        0        0   1.0000
          0   1.0000  -1.0000        0        0   1.0000
       
When combining the real poles and zeros we get
  sos2 =
     7.9956e-04   1.5991e-03   7.9956e-04   1.0000e+00  -1.9263e+00   9.2963e-01
     1.0000e+00  -2.0000e+00   1.0000e+00   1.0000e+00  -1.9914e+00   9.9140e-01       


Mathematically, both solutions seem to be correct, because when recalculating the transfer functions, we get same result

  >> [B,A]=sos2tf(sos)
  B =
     7.9956e-04            0  -1.5991e-03            0   7.9956e-04
  A =
     1.0000  -3.9176   5.7569  -3.7609   0.9216
 
  >> [B,A]=sos2tf(sos2)
  B =
     7.9956e-04            0  -1.5991e-03            0   7.9956e-04
  A =
     1.0000  -3.9176   5.7569  -3.7609   0.9216


The problem is the function sosfilt, which seems not to handle transfer functions with leading zeros in the denominator.

To overcome your problem:
- edit the function  \octave-7.3.0-w64\mingw64\share\octave\packages\signal-1.4.3\zp2sos.m
- change the value of DoNotCombineReal to 0 on line 78.

For the future, I will reexamine the choice of this parameter:
Value to be modified?
Value depending on filter type?
Default value rather that fixed value?
etc.

Suggestions are welcome.


Charles Praplan <charprap>
Tue 29 Nov 2022 12:18:17 PM UTC, original submission:  


sampling_rate = 44100;
length = 30;
hpf = 45;
lpf = 450;
filter_order = 2;
typenoise = noise((length/2)*sampling_rate, 1, 'pink');
[z, p, k] = butter(filter_order, [hpf/(sampling_rate/2), lpf/(sampling_rate/2)]);
sos = zp2sos (z, p, k);
filtered = sosfilt(sos, typenoise);


"filtered" is resulting a NaN

Same code as working fine on Octave 6.4

Renato <yamane>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #54190:  slope.jpg added by yamane (182KiB - image/jpeg)
file #54180:  Histogram_original.png added by yamane (23KiB - image/png)
file #54164:  zp2sos.m added by charprap (6KiB - text/plain)
file #54165:  sos2tf.m added by charprap (5KiB - text/plain)
file #54166:  zp2sos.diff added by charprap (1KiB - application/octet-stream)
file #54167:  sos2tf.diff added by charprap (744B - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by lostbard (Posted a comment)
  • -email is unavailable- added by lostbard (adding signal maintainer )
  • -email is unavailable- added by charprap (Posted a comment)
  • -email is unavailable- added by yamane (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 17 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-05-17 lostbard StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2023-04-25 lostbard Item GroupUnexpected Error or Warning Incorrect Result
        StatusNone Ready For Test
        Assigned toNone lostbard
        Release7.3.0 other
        Operating SystemMicrosoft Windows Any
    2023-04-21 lostbard Carbon-Copy- Added mtmiller
    2023-01-04 yamane Attached File- Added slope.jpg, #54190
    2023-01-03 yamane Attached File- Added Histogram_original.png, #54180
        Attached File- Added Histogram_proposed_change.png, #54181
    2022-12-28 charprap Attached File- Added zp2sos.m, #54164
        Attached File- Added sos2tf.m, #54165
        Attached File- Added zp2sos.diff, #54166
        Attached File- Added sos2tf.diff, #54167
    2022-11-29 mmuetzel Carbon-CopyRemoved 102357 -
    2022-11-29 mmuetzel Summarysosfilt - NaN [octave forge] (signal) sosfilt - NaN

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code