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

 
 

bug #53855: [octave forge] (signal) document that IIR filters in transfer function form may suffer round-off errors

Submitter:  Renato <yamane>
Submitted:  Mon 07 May 2018 10:07:04 AM UTC
   
 
Category:  Octave Package Severity:  2 - Minor
Priority:  3 - Low Item Group:  Documentation
Status:  Postponed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 08 May 2018 08:00:41 PM UTC, comment #7: 

Q: This bug remains with the question: is any description needed in the butter (and cheby1, cheby2, ellip) doc string to mention this problem?

A: Of course, and an workaround will be appreciated :-)

Renato <yamane>
Tue 08 May 2018 07:45:19 PM UTC, comment #6: 

Those are questions for the help mailing list, not the bug tracker.

The freqz function in Octave does not yet support a filter specified by sos (second-order sections), and the fvtool is not yet implemented either. The only way to visualize the filter frequency response at the moment is with freqz, which only supports transfer function format. A wish list bug requesting support for second-order sections in the freqz function would be appropriate.

The filter function only supports transfer function format, you would use the sosfilt function to apply a second-order section filter to your 'typenoise' signal.

This bug remains with the question: is any description needed in the butter (and cheby1, cheby2, ellip) doc string to mention this problem?

Mike Miller <mtmiller>
Group Member
Tue 08 May 2018 06:59:28 PM UTC, comment #5: 

The MatLab documentation is not very well intuitive :-(

The documentation also says to use "hfvt", but it´s not implemented in Octave.

How to use ZPK Design instead a TF Design of this code below? :-)


n = 5;
hpf = 50;
lpf = 10*hpf;
sampling_rate = 44100;

typenoise = noise(60*sampling_rate, 1, 'pink');
[b,a] = butter(n, [hpf/(sampling_rate/2), lpf/(sampling_rate/2)]);
filtered = filter(b, a, typenoise);

freqz (b, a, 4096, sampling_rate)
ax = findall (gcf, 'type', 'axes');
set (ax, 'xlim', [hpf/4 lpf*4]);
set (ax, 'xscale', 'log');
subplot (2, 1, 1);
set (gca, 'ylim', [-35 3]);


Thank you,
Renato

Renato <yamane>
Tue 08 May 2018 06:02:58 PM UTC, comment #4: 

Agree, this is exactly shown in the example on the Matlab help page.

The actual computed coefficients are consistent with both Matlab and SciPy.

Do we need to document this in the Octave butter function? Any doc string patches or suggestions?

Mike Miller <mtmiller>
Group Member
Tue 08 May 2018 04:04:13 PM UTC, comment #3: 

Is not this the numerical instability as described in:


https://www.mathworks.com/help/signal/ref/butter.html

<<<
In general, use the [z,p,k] syntax to design IIR filters. To analyze or implement your filter, you can then use the [z,p,k] output with zp2sos

>>>


Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 08 May 2018 03:18:23 PM UTC, comment #2: 

Same problem in MatLab...
File attached.

IF you increase the high frequency cut-off, the result is much better. But I didn't undestand why...

In the example below, I changed the high-frequency cut-off from 100Hz to 200Hz.


n = 5;
hpf = 200;
lpf = 10*hpf;
[b, a] = butter (n, [hpf/22050, lpf/22050]);
freqz (b, a, 4096, 44100)
ax = findall (gcf, 'type', 'axes');
set (ax, 'xlim', [hpf/4 lpf*4]);
set (ax, 'xscale', 'log');
subplot (2, 1, 1);
set (gca, 'ylim', [-35 3]);



Renato <yamane>
Mon 07 May 2018 06:53:23 PM UTC, comment #1: 

Confirmed here. It would be useful to compare the filter coefficients to other software (e.g. SciPy).

Here is a simpler example, since I'm not sure where the 'noise' and 'normalize' functions are coming from, and they aren't needed.


n = 5;
[b, a] = butter (n, [100/22050, 1000/22050]);
freqz (b, a, 4096, 44100);
ax = findall (gcf, 'type', 'axes');
set (ax, 'xlim', [20 4000]);
set (ax, 'xscale', 'log');
subplot (2, 1, 1);
set (gca, 'ylim', [-35 3]);


Mike Miller <mtmiller>
Group Member
Mon 07 May 2018 10:07:04 AM UTC, original submission:  


%=== example of the code ===
n=5;
typenoise = noise(60*44100, 1, 'pink');

[b,a] = butter(n, [100/(22050), 1000/(22050)]);
filtered = filter(b, a, typenoise);
normalized = normalize(filtered, 'wav');
audiowrite ('test.wav', normalized, 44100);

freqz(b, a, 4096, 44100);
ax = findall(gcf, 'Type', 'axes');
set(ax, 'XLim', [20 4000]);
set(ax, 'YLim', [-35 3]);
set(ax, 'XScale', 'log');
%=============================


When n<=4, the result in low-frequency is OK.
When n >=5, the result in low-frequency is strange.

Renato <yamane>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #44120:  test_RYamane.png added by yamane (55KiB - image/png - MatLab r2017b)
file #44109:  n5.png added by yamane (10KiB - image/png)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by dasergatskov (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 13 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-27 mtmiller Carbon-CopyRemoved 80942 -
    2018-06-05 mtmiller Item GroupInaccurate Result Documentation
        StatusConfirmed Postponed
        Summary[octave forge] (signal) butter band pass filter with higher order has round-off errors [octave forge] (signal) document that IIR filters in transfer function form may suffer round-off errors
    2018-05-08 mtmiller Severity3 - Normal 2 - Minor
    2018-05-08 mtmiller Priority5 - Normal 3 - Low
        Summary[octave forge] (signal) butter band pass filter with higher order has strange pass band [octave forge] (signal) butter band pass filter with higher order has round-off errors
    2018-05-08 yamane Attached File- Added test_RYamane.png, #44120
    2018-05-07 mtmiller Item GroupIncorrect Result Inaccurate Result
        StatusNone Confirmed
        Release4.4.0 dev
        Summary[signal] butter filter with strange results [octave forge] (signal) butter band pass filter with higher order has strange pass band
    2018-05-07 yamane Attached File- Added n5.png, #44109

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code