# Sample rate severely influences quality of ButterworthFilterModel?

I've hesitated for a while about whether this should be posted to dsp.SE or not, but still decide to discuss it here.

Today I was playing with ButterworthFilterModel in order to answer this question but found something strange, that is, the filtering quality of the same ButterworthFilterModel seems to be influenced severely by the sample rate. Consider the following example, whose target is filtering out component with frequency higher than $15\text{Hz}$ from the signal $\sin(20\pi t)+\sin(40\pi t)$:

help[sr_, modelfunc_: ButterworthFilterModel] := (
tend = 10;
time = Range[0, tend, 1/sr];

freq@1 = 10; freq@2 = 20;
data = Sin[2 Pi freq@1 time] + Sin[2 Pi freq@2 time];
fp = 12; fs = 14; ap = 1.; as = 20.;
model = modelfunc[{"Lowpass", {fp, fs} 2 Pi, {ap, as}}];
dismodel = ToDiscreteTimeModel[model, 1/sr];
dataafter = RecurrenceFilter[dismodel, data]
)

ListLinePlot[help[10^2], DataRange -> {0, tend}, PlotRange -> {{1, 2}, Automatic}]


As we can see, the filter works quite well with sample rate $100$, but, if I modify the sample rate to $1000$:

help[10^3] // Max
(* 1.228003891092094*10^1083 *)


The filter no longer works properly. Furthermore, if I turn to another filter model e.g.:

ListLinePlot[help[#, Chebyshev1FilterModel] & /@ {10^2, 10^3}, DataRange -> {0, tend},
PlotRange -> {{1, 2}, Automatic}]


The filter works properly now.

I've also tested in MATLAB/Octave, the result is similar:

function [time,dataafter]=help(sr)
tend=10;
time=0:1/sr:tend;
freq1=10;freq2=20;
data=sin(2*pi*freq1*time)+sin(2*pi*freq2*time);
fp=12;fs=14;ap=1;as=20;
[n,Wn]=buttord(fp*2/sr,fs*2/sr,ap,as);
[b,a]=butter(n,Wn);
dataafter=filter(b,a,data);
end

[x,dat]=help(100);
plot(x(100:200),dat(100:200))


[x,dat]=help(1000);
max(dat)
# 4.0522e+301


So this seems to be the nature of Butterworth filter model.

My question is:

1. Why does this happen?

2. If the first question is too hard, is there a systematic way to decide proper parameters for Butterworth filter?

  sr = 1000.;
tend = 10.;
time = Range[0, tend, 1/sr];
freq@1 = 10.; freq@2 = 20.;
data = Sin[2 Pi freq@1 time] + Sin[2 Pi freq@2 time];
ListLinePlot[Transpose[{time, data}][[1 ;; 300]]]


Data looks fine. Note I have made everything reals. I don't want calculations done in integers. I assume you are trying to get rid of the 20 Hz sine wave.

fp = 12.; fs = 14.; ap = 1.; as = 20.;
model = ButterworthFilterModel[{"Lowpass", {fp, fs} 2 Pi, {ap, as}}];
dismodel = ToDiscreteTimeModel[model, 1/sr];
dataafter = RecurrenceFilter[dismodel, data];
MinMax@dataafter


{-3.447895715878780*10^1078, 5.294173234545417*10^1083}

Big numbers. So what has happened here is that your filter has become unstable. This happens when you have very demanding specifications.

Let's try a notch filter.

fp = 19.5; fs = 20.5;
model = ButterworthFilterModel[{"BandStop", 3, {fp, fs} 2 Pi}];
dismodel = ToDiscreteTimeModel[model, 1/sr];
dataafter = RecurrenceFilter[dismodel, data];

n1 = 1; n2 = 1200;
ListLinePlot[{data[[n1 ;; n2]], Sin[2 Pi freq@1 time[[n1 ;; n2]]],
dataafter[[n1 ;; n2]]}]


The filter output is taking awhile to settle down. This happens when the bandwidth is small. Let's look at a later time.

n1 = 1200; n2 = 2400;
ListLinePlot[{data[[n1 ;; n2]], Sin[2 Pi freq@1 time[[n1 ;; n2]]],
dataafter[[n1 ;; n2]]},
PlotLegends -> LineLegend[{"Data", "10 Hz", "Filtered Data"}]]


This looks reasonable. Hope that helps.

Edit

The issue you are seeing is a result of polynomials with large coefficients. It happens like this. Filters are based on analogue devices that have an input x(t) and an output y(t). A typical relationship between the input and output is

$c_0 y+c_1 \dot{y}+\ddot{y}=b_2 \ddot{x}+b_0 x+b_1 \dot{x}$

If this is discretized in the usual way by making

$\dot{y}=\frac{y_{n+1}-y_n}{T}$

where T is the time step then you get a difference equation of the form

$y_n=\alpha _0 x_n+\alpha _1 x_{n-1}+\alpha _2 x_{n-2}-\beta _1 y_{n-1}-\beta _2 y_{n-2}$

where the coefficients α and β all have terms with T or $T^2$ on the denominator. The time step is the reciprocal of the sample rate so that for large sample rates T is very small and the coefficients are very large or very small. Thus the difference equation above has to deal with adding terms with big and small coefficients and the usual problem of loss of numerical accuracy occurs. Note also that there is feedback with previous output values interacting with the current output thus allowing for instability.

Here is a module that extracts the coefficients from a Butterworth filter.

    ClearAll[findCoefficients];
findCoefficients[f1_, f2_, ap_, as_, sr_] :=
Module[{bpf, dtm, cnum, cden},
bpf = ButterworthFilterModel[{"Lowpass",
2 \[Pi] {f1, f2}, {ap, as}}];
dtm = ToDiscreteTimeModel[bpf, 1/sr];
cnum = Reverse@CoefficientList[dtm[[1, 1, 1, 1]], \[FormalZ]];
cden = Reverse@CoefficientList[dtm[[1, 2, 1, 1]], \[FormalZ]];
{cnum, cden}
]


If we apply this to your first case we get

f1 = 12.; f2 = 14.; ap = 1; as = 20.;
{cnum, cden} = Re@findCoefficients[f1, f2, ap, as, 100]


{{0.00692843, 0.138569, 1.3164, 7.89841, 33.5683, 107.418, 268.546, 537.092, 872.775, 1163.7, 1280.07, 1163.7, 872.775, 537.092, 268.546, 107.418, 33.5683, 7.89841, 1.3164, 0.138569, 0.00692843},

{1.39374*10^8, -1.46698*10^9, 7.52801*10^9, -2.49629*10^10, 5.98438*10^10, -1.10034*10^11, 1.6075*10^11, -1.90814*10^11, 1.86697*10^11, -1.51903*10^11, 1.0325*10^11, -5.86861*10^10, 2.78257*10^10, -1.09393*10^10, 3.52908*10^9, -9.19423*10^8, 1.88822*10^8, -2.94484*10^7, 3.27984*10^6, -232522., 7888.92}}

Here the first coefficients are for the α terms and the second set of coefficients are for the β terms. There are 21 coefficients for each implying a differential equation of order 20. The coefficients are reasonable although the β ones are large with alternating signs. This will lead to only a few digits of accuracy remaining after each sum in the recurrence equation has been completed.

Now look at what happens if we increase the sample rate.

f1 = 12.; f2 = 14.; ap = 1; as = 20.;
{cnum, cden} = Re@findCoefficients[f1, f2, ap, as, 1000]


{{6.92843*10^-23, 1.38569*10^-21, 1.3164*10^-20, 7.89841*10^-20, 3.35683*10^-19, 1.07418*10^-18, 2.68546*10^-18, 5.37092*10^-18, 8.72775*10^-18, 1.1637*10^-17, 1.28007*10^-17, 1.1637*10^-17, 8.72775*10^-18, 5.37092*10^-18, 2.68546*10^-18, 1.07418*10^-18, 3.35683*10^-19, 7.89841*10^-20, 1.3164*10^-20, 1.38569*10^-21, 6.92843*10^-23},

{1.72349*10^6, -3.27575*10^7, 2.95778*10^8, -1.68698*10^9, 6.8163*10^9, -2.074*10^10, 4.93078*10^10, -9.37932*10^10, 1.4498*10^11, -1.83902*10^11, 1.92476*10^11, -1.66509*10^11, 1.18852*10^11, -6.96177*10^10, 3.31367*10^10, -1.26196*10^10, 3.75513*10^9, -8.41435*10^8, 1.3357*10^8, -1.33929*10^7, 637957.}}

Now the α terms are very small. If Chop was used they would all be zero. This will not give good answers.

We can now look at what happens to the coefficients as the sample rate is incresaed

f1 = 12.; f2 = 14.; ap = 1; as = 20.;
coffs = Table[{sr, #} & /@
Abs@Flatten@findCoefficients[f1, f2, ap, as, sr], {sr,
Range[10] 100}];
ListLogPlot[coffs, PlotRange -> All, Frame -> True,
FrameLabel -> {"Sample Rate", "Coefficients"}]


As you can see the coefficients have very large and very small values so that the usual problem of loss of precision occurs. I did look at using higher precision here with a RecurrenceFilter and it does help.

• Values of ap and as are not important, we can adjust them freely, the only requirement is the $20\text{Hz}$ signal should be filter out by low pass filter. Well, actually obtaining such a filter is not too hard by trial and error. As mentioned in the question, I'm mainly interested in, can this be done systematically? Say, we have a signal, and I managed to figure out a set of proper parameters for the filter under a certain sample rate, can we deduce the proper parameters for another sample rate then? – xzczd Jun 28 '18 at 11:38