Filter Offsets Gains?

[mcbaron]

I'm making a filterbank to do some signal analysis. I've constructed an array of bandpass filters, which reconstruct. Here are the frequency responses of three adjacent bins:

If I take some real-world data, and run it through those same 3 filters, I get a gain offset between each of them.

y1 is constructed: X = periodogram(x, fs=fs) plot(freq(X), pow2db.(power(X))) and y2 through y4 are constructed using a = filt(fb[4], x) A = periodogram(a, fs=fs) plot!(freq(A), pow2db.(power(A)))

I'm really confused what is going wrong here. I'd expect the filtered spectrums to have all roughly the same gain as eachother, with the passband to be about the same level as the original signal.

[galenlynch]

You didn't post enough information to determine what happened in your case, but I can't replicate your issue:

using DSP, PyPlot
a = rand(1000) .- 0.5;
b1 = digitalfilter(Bandpass(0.1, 0.2), FIRWindow(blackman(151)));
b2 = digitalfilter(Bandpass(0.4, 0.5), FIRWindow(blackman(151)));
c1 = filtfilt(b1, a);
c2 = filtfilt(b2, a);
s_raw = periodogram(a);
s_filt1 = periodogram(c1);
s_filt2 = periodogram(c2);
plot(freq(s_raw), pow2db.(power(s_raw)), label = "input");
plot(freq(s_filt1), pow2db.(power(s_filt1)), label = "first_bandpass");
plot(freq(s_filt2), pow2db.(power(s_filt2)), label = "second_bandpass");
legend();

[mcbaron]

Generated with the following code:

# using linkwitz-reily (squared butterworth)
using DSP
x = x1[:, 36, 1]
filtermethod = Butterworth(20)
hp1 = digitalfilter(Highpass(180, fs=44100), filtermethod)
hp1 = hp1 * hp1
lp1 = digitalfilter(Lowpass(226, fs=44100), filtermethod)
lp1 = lp1 * lp1
hp2 = digitalfilter(Highpass(226, fs=44100), filtermethod)
hp2 = hp2 * hp2
lp2 = digitalfilter(Lowpass(285, fs=44100), filtermethod)
lp2 = lp2 * lp2
b1 = hp1 * lp1
b2 = hp2 * lp2
c1 = filt(b1, x)
c2 = filt(b2, x)
X = periodogram(x, fs=44100)
C1 = periodogram(c1, fs=44100)
C2 = periodogram(c2, fs=44100)
plot(freq(X), pow2db.(power(X)), xscale=:log10, xlims=(10, 20000), label="input")
plot!(freq(C1), pow2db.(power(C1)), label="firstpass")
plot!(freq(C2), pow2db.(power(C2)), label="secondpass")

[galenlynch]

Could you please verify that problem is with DSP.jl, and not your filter design, by seeing if other DSP packages such as SciPy produce different results? I don't know enough about Linkwitz-Reily filters to tell exactly what is happening with your example, but it seems a little unusual and I don't know what results I'd expect. I'm only familiar with two-way crossover designs that use a high-pass and low-pass filter, each constructed with a cascade of two first order Butterworth filters. In your example, you are making a multi-way design with band-pass filters constructed by combining low-pass and high-pass filter, each constructed with a cascade of 20-order Butterworth filters. I'm not seeing any problem if you use band-pass filters with more straight-forward designs:

julia> using DSP, PyPlot

julia> a = rand(10000) .- 0.5;

julia> function show_periodogram(signal, signal_name)
       psd = periodogram(signal)
       plot(freq(psd), pow2db.(power(psd)), label = signal_name)
       end
show_periodogram (generic function with 1 method)

julia> bp1 = digitalfilter(Bandpass(0.2, 0.3), Butterworth(2));

julia> bp2 = digitalfilter(Bandpass(0.6, 0.7), Butterworth(2));

julia> show_periodogram(a, "input")
Warning: Ignoring XDG_SESSION_TYPE=wayland on Gnome. Use QT_QPA_PLATFORM=wayland to run on Wayland anyway.
1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7f6b1c7caba8>

julia> show_periodogram(filt(bp1, a), "first band-pass")
1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7f6b1c665d68>

julia> show_periodogram(filt(bp2, a), "second band-pass")
1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7f6b1c64ab70>

julia> legend()
PyObject <matplotlib.legend.Legend object at 0x7f6b1c683e10>

[galenlynch]

Yeah when I try your filters (constructed by copying your example above) on just white noise (x = rand(100_000) .- 0.5;) it looks they don't have unity gain:

[mcbaron]

In my first post I graphed the frequency response of each of those filters, plot(f, amp2db.(abs.(freqz(filter,f,fs=44100)))) Which shows they all have unity gain.

(I'm posting from my phone, so apologies for incorrect code formatting)

I also have some code that sums all the filters in the bank, which reconstructs perfectly to a unity gain band pass, at least according to the graph. I'll post code and plots when I'm back at my desk.

On Sat, Jun 15, 2019, 4:59 PM Galen Lynch notifications@github.com wrote:

Yeah when I try your filters (constructed by copying your example above) on just white noise (x = rand(100_000) .- 0.5;) it looks they don't have unity gain: [image: Figure_1] https://user-images.githubusercontent.com/6831311/59556352-ec0e5100-8f8e-11e9-9465-3a1bbb7b47bb.png

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaDSP/DSP.jl/issues/309?email_source=notifications&email_token=ABAGHMVNN3FFLZ3MVE7PNJTP2VJ4PA5CNFSM4HYL2GT2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODXZAAXY#issuecomment-502399071, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAGHMS3N5OLJPWRHVQMJZDP2VJ4PANCNFSM4HYL2GTQ .

[galenlynch]

I suppose you're right with the frequency response. If you look closely at the frequency response of your filters it seems like the pass band may still be somewhat attenuated, and are very narrow for both filters.

julia> using DSP, PyPlot

julia> function show_frequency_response(b, frs, fs, label)
       r = freqz(b, frs, fs)
       semilogx(frs, amp2db.(abs.(r)), label = label)
       end

julia> function show_periodogram(signal, fs, signal_name)
              psd = periodogram(signal,fs = fs)
              semilogx(freq(psd), pow2db.(power(psd)), label = signal_name)
              end

julia> filtermethod = Butterworth(20);

julia> hp1 = digitalfilter(Highpass(180, fs=44100), filtermethod);

julia> hp1 = hp1 * hp1;

julia> lp1 = digitalfilter(Lowpass(226, fs=44100), filtermethod);

julia> lp1 = lp1 * lp1;

julia> hp2 = digitalfilter(Highpass(226, fs=44100), filtermethod);

julia> hp2 = hp2 * hp2;

julia> lp2 = digitalfilter(Lowpass(285, fs=44100), filtermethod);

julia> lp2 = lp2 * lp2;

julia> b1 = hp1 * lp1;

julia> b2 = hp2 * lp2;

julia> show_frequency_response(b1, 0:1:1000, 44100, "first band-pass");

julia> show_frequency_response(b2, 0:1:1000, 44100, "second band-pass"); # Zoomed-out filter response

julia> legend();

julia> xlim(184, 280);

julia> ylim(-1, 0);  # zoomed in filter response

julia> figure();

julia> show_periodogram(a, 44100, "input");

julia> show_periodogram(filt(b1, a), 44100, "first band-pass");

julia> show_periodogram(filt(b2, a), 44100, "second band-pass"); 

julia> legend(); # filter results

I can't even re-create this example in Matlab, possibly because the transfer function has both enormous and miniscule coefficients. Here was my attempt:

fs = 44100;
[zlp,plp,klp] = butter(20, 2 * 286 / fs,'low');
[zhp,php,khp] = butter(20, 2 * 226 / fs,'high');
zcasc = [zlp; zlp; zhp; zhp];
pcasc = [plp; plp; php; php];
kcasc = klp ^2 * khp ^2;
[b,a] = zp2tf(zcasc, pcasc, kcasc);
nsig = 10000;
x = rand(nsig, 1) - 0.5;
win = kaiser(nsig, 38);
[praw, fraw] = periodogram(x, win, fs);
semilogx(fraw, pow2db(praw), 'DisplayName', 'input');
semilogx(fraw, pow2db(periodogram(filter(b, a, x), win, fs)), 'DisplayName', 'band-passed')

If you look at the output of filter(b,a,x) then you see it's just NaNs.

Looking at the PolynomialRatio form of your filters in Julia shows that this might be the underlying problem:

julia> PolynomialRatio(b1)
PolynomialRatio{Float64}(Poly(8.996718225380656e-73 - 3.5986872901522623e-71*x^2 + 7.017440215796912e-70*x^4 - 8.888757606676088e-69*x^6 + 8.222100786175382e-68*x^8 - 5.919912566046275e-67*x^10 + 3.45328233019366e-66*x^12 - 1.6773085603797778e-65*x^14 + 6.918897811566583e-65*x^16 - 2.460052555223674e-64*x^18 + 7.626162921193389e-64*x^20 - 2.0798626148709243e-63*x^22 + 5.026334652604735e-63*x^24 - 1.0825951559456351e-62*x^26 + 2.087862086466582e-62*x^28 - 3.6189609498754086e-62*x^30 + 5.654626484180326e-62*x^32 - 7.983002095313402e-62*x^34 + 1.0200502677344901e-61*x^36 - 1.1811108363241465e-61*x^38 + 1.2401663781403539e-61*x^40 - 1.1811108363241465e-61*x^42 + 1.0200502677344901e-61*x^44 - 7.983002095313402e-62*x^46 + 5.654626484180326e-62*x^48 - 3.6189609498754086e-62*x^50 + 2.087862086466582e-62*x^52 - 1.0825951559456351e-62*x^54 + 5.026334652604735e-63*x^56 - 2.0798626148709243e-63*x^58 + 7.626162921193389e-64*x^60 - 2.460052555223674e-64*x^62 + 6.918897811566583e-65*x^64 - 1.6773085603797778e-65*x^66 + 3.45328233019366e-66*x^68 - 5.919912566046275e-67*x^70 + 8.222100786175382e-68*x^72 - 8.888757606676088e-69*x^74 + 7.017440215796912e-70*x^76 - 3.5986872901522623e-71*x^78 + 8.996718225380656e-73*x^80), Poly(0.2288697481761285 - 18.64705527007814*x + 750.1378202550235*x^2 - 19863.233534845487*x^3 + 389419.82042536384*x^4 - 6.02838307409858e6*x^5 + 7.674531791938004e7*x^6 - 8.262822775345889e8*x^7 + 7.679023819184682e9*x^8 - 6.256659800573764e10*x^9 + 4.524277830199429e11*x^10 - 2.932275932736003e12*x^11 + 1.7172152388154312e13*x^12 - 9.148379822907816e13*x^13 + 4.459094624102793e14*x^14 - 1.9982827559423302e15*x^15 + 8.268176392811268e15*x^16 - 3.170312319054675e16*x^17 + 1.1301424805149573e17*x^18 - 3.7560921737015386e17*x^19 + 1.166817355283205e18*x^20 - 3.395498033076769e18*x^21 + 9.274773449862842e18*x^22 - 2.3821867979232027e19*x^23 + 5.762535632943463e19*x^24 - 1.3147333827987443e20*x^25 + 2.8327283743099473e20*x^26 - 5.770516443107177e20*x^27 + 1.1125361467887413e21*x^28 - 2.0319028563933913e21*x^29 + 3.5183349597356126e21*x^30 - 5.780064347362568e21*x^31 + 9.015043211432437e21*x^32 - 1.3356295150842e22*x^33 + 1.880605520178796e22*x^34 - 2.517573239337535e22*x^35 + 3.205447062040952e22*x^36 - 3.882736401006067e22*x^37 + 4.475308485795594e22*x^38 - 4.909188820154042e22*x^39 + 5.125515424905284e22*x^40 - 5.093537290560219e22*x^41 + 4.8177303190096e22*x^42 - 4.3367778479899815e22*x^43 + 3.7147333358360558e22*x^44 - 3.0271289500221416e22*x^45 + 2.3461528260198684e22*x^46 - 1.728838379287955e22*x^47 + 1.210726393750307e22*x^48 - 8.054168230873824e21*x^49 + 5.08668742466375e21*x^50 - 3.047969891463747e21*x^51 + 1.7315368899195579e21*x^52 - 9.318418510960093e20*x^53 + 4.7461602563845266e20*x^54 - 2.2855204991807704e20*x^55 + 1.039372087432625e20*x^56 - 4.458032872227673e19*x^57 + 1.8008636200149017e19*x^58 - 6.840551167264113e18*x^59 + 2.438937261911148e18*x^60 - 8.145995508202497e17*x^61 + 2.543029179953258e17*x^62 - 7.40168255424158e16*x^63 + 2.0028501131560204e16*x^64 - 5.022338893577209e15*x^65 + 1.1628030709180018e15*x^66 - 2.475221845948516e14*x^67 + 4.820645444481492e13*x^68 - 8.540747130022041e12*x^69 + 1.3672594191279043e12*x^70 - 1.9618007156927527e11*x^71 + 2.4982109892179073e10*x^72 - 2.7890888818568125e9*x^73 + 2.687798389254786e8*x^74 - 2.1905672647637635e7*x^75 + 1.4681986690841657e6*x^76 - 77701.2163443763*x^77 + 3044.5989981323955*x^78 - 78.52546953593814*x^79 + 1.0*x^80))

So you might be running into numerical instability.

At any rate, I think the problem is with your filters, and not DSP.jl. Matlab gives up before DSP.jl.

[galenlynch]

The issue is not just one of filter order, however. If you make a band-pass filter of the same order and passband as in your example with a simple bilinear transform, the filtering results look quite reasonable:

using DSP, PyPlot
function show_periodogram(signal, fs, signal_name)
              psd = periodogram(signal,fs = fs)
              semilogx(freq(psd), pow2db.(power(psd)), label = signal_name)
              end
a = rand(100000) .- 0.5;
bpzpk = digitalfilter(Bandpass(226, 285, fs = 44100), Butterworth(40));
show_periodogram(a, 44100, "input");
show_periodogram(filt(bpzpk, a), 44100, "bp");
legend();

[mcbaron]

Thank you for your analysis!

It raises two questions in my mind however: 1) I recreated your example with 40th order bandpass filters, and am still seeing some level offset between the two filtered responses.

The resultant filters don't reconstruct exactly as desired:

graph_freq_resp(b1, 10:22000, fs, "band1")
graph_freq_resp!(b2, 10:22000, fs, "band2")
plot!(10:22000, amp2db.(abs.(freqz(b1, 10:22000, fs) + freqz(b2, 10:22000, fs))), label="sum")

Avoiding the 3db bump at the crossover is the reason for using squared butterworth filters.

2) Jumping back to the way I was making Linkwitz-Reily filters, since there is an offset, is there a way to normalize the resultant filters to be unity gain, possibly in the process zeroing out the very small filter coefficients?

I've noticed that the fields of a filter object are not directly mutable ie I can't do an assignment like filter.k = 2*filter.k to double the gain of a ZPK object. If I wanted to take my designed filter to a PolinomialRatio then threshold the coefficients, how would I do that?

[galenlynch]

Not sure what's happening in your case with your data, both BP filters seem fine with white noise input:

using DSP, PyPlot
function show_periodogram(signal, fs, signal_name)
              psd = periodogram(signal,fs = fs)
              semilogx(freq(psd), pow2db.(power(psd)), label = signal_name)
              end
a = rand(100000) .- 0.5;
bpzpk1 = digitalfilter(Bandpass(180, 226, fs = 44100), Butterworth(40));
bpzpk2 = digitalfilter(Bandpass(226, 285, fs = 44100), Butterworth(40));
show_periodogram(a, 44100, "input");
show_periodogram(filt(bpzpk1, a), 44100, "bp1");
show_periodogram(filt(bpzpk2, a), 44100, "bp2");
legend();

You can directly construct a ZeroPoleGain or PolynomialRatio with arbitrary poles etc. as described in the documentation. You can also convert between them by calling their constructors.

I just verified that you can construct more "traditional" Linkwitz-Reily filters (e.g. a second order two-way cross-over) and not run into these numerical stability problems, so maybe you just need to back off on the order of your filters.

I'm closing this issue because it's not a problem with the package, but instead with its application. I suspect that your issues stem from numerical instability. I think you might just need to decrease your filter order. I would recommend posting on Discourse to get more help with your questions. Feel free to re-open it if you can isolate a similar problem with a minimal working example, and ideally can verify that your minimal working example works as intended in other DSP packages.

[mcbaron]

Thank you again, this has been a really helpful discussion for me.

[galenlynch]

My pleasure! If you do post this on Discourse I might be able to chime in over there as well.