sosfit in filter function

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

sosfit in filter function

yamane
Hi,

There is a numerical instability on butter function, as described in
https://www.mathworks.com/help/signal/ref/butter.html (topic "Limitations").

Is recommend use [z, p, k] syntax and sosfilt function to apply a
second-order section filter, as described in
https://savannah.gnu.org/bugs/?53855

But, the documentation of SOSFILT is not clear:
https://octave.sourceforge.io/signal/function/sosfilt.html

Someone have an idea about how to use ZPK Design instead a TF Design in 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);
audiowrite (filename.wav, filtered, sampling_rate);

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]);
-----------------------------

Best regards,
Renato



--
Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sosfit in filter function

Mike Miller-4
On Tue, May 08, 2018 at 23:21:33 -0700, yamane wrote:

> There is a numerical instability on butter function, as described in
> https://www.mathworks.com/help/signal/ref/butter.html (topic "Limitations").
>
> Is recommend use [z, p, k] syntax and sosfilt function to apply a
> second-order section filter, as described in
> https://savannah.gnu.org/bugs/?53855
>
> But, the documentation of SOSFILT is not clear:
> https://octave.sourceforge.io/signal/function/sosfilt.html
>
> Someone have an idea about how to use ZPK Design instead a TF Design in this
> code below?
Instead of this

    [b, a] = butter (...
    output = filter (b, a, input);

you would use this

    [z, p, k] = butter (...
    sos = zp2sos (z, p, k);
    output = sosfilt (sos, input);

--
mike



signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: sosfit in filter function

yamane
2018-05-09 19:31 GMT+02:00 Mike Miller <[hidden email]>:

> Instead of this
>
>     [b, a] = butter (...
>     output = filter (b, a, input);
>
> you would use this
>
>     [z, p, k] = butter (...
>     sos = zp2sos (z, p, k);
>     output = sosfilt (sos, input);

Thank you so much for your support!
Now, as it´s not possible see the response frequency of the filter
using freqz, I opened the output file in Audacity.

When the audio is very compressed (6dB Crest Factor), the response
frequency don´t keep attenuating in LOW frequencies.
See here: https://ibb.co/iQD8Ty

++++++++++++++++++++++++++++
sampling_rate = 44100;
length = 120;
hpf = 40;
lpf = 400;
filter_order = 4;
crest_factor = 6;

typenoise = noise((length)*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 (normalized(normalized > 1) || normalized(normalized < -1))
  normalized(normalized > 1) = 1;
  normalized(normalized < -1) = -1;
  normalized = normalized / (rms(normalized) / 10^(-crest_factor/20));
endwhile

audiowrite ('output.wav', normalized, sampling_rate);
----------------------------------------------------

Best regards,
Renato


Reply | Threaded
Open this post in threaded view
|

Re: sosfit in filter function

yamane
2018-05-09 20:19 GMT+02:00 Renato S. Yamane <[hidden email]>:

> 2018-05-09 19:31 GMT+02:00 Mike Miller <[hidden email]>:
>> Instead of this
>>
>>     [b, a] = butter (...
>>     output = filter (b, a, input);
>>
>> you would use this
>>
>>     [z, p, k] = butter (...
>>     sos = zp2sos (z, p, k);
>>     output = sosfilt (sos, input);
>
> Thank you so much for your support!
> Now, as it´s not possible see the response frequency of the filter
> using freqz, I opened the output file in Audacity.
>
> When the audio is very compressed (6dB Crest Factor), the response
> frequency don´t keep attenuating in LOW frequencies.
> See here: https://ibb.co/iQD8Ty
>
> ++++++++++++++++++++++++++++
> sampling_rate = 44100;
> length = 120;
> hpf = 40;
> lpf = 400;
> filter_order = 4;
> crest_factor = 6;
>
> typenoise = noise((length)*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 (normalized(normalized > 1) || normalized(normalized < -1))
>   normalized(normalized > 1) = 1;
>   normalized(normalized < -1) = -1;
>   normalized = normalized / (rms(normalized) / 10^(-crest_factor/20));
> endwhile
>
> audiowrite ('output.wav', normalized, sampling_rate);
> ----------------------------------------------------

As an workaround, I tried use pwelch:
pwelch(normalized,hanning(8192), 'loglog');

The "shape" of the graphic looks OK, but the axis are a little bit strangers.
Can someone give me a help, again? :-)

Thank you,
Renato