8
$\begingroup$

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}]

Mathematica graphics

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}]

Mathematica graphics

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

enter image description here

[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?

$\endgroup$
9
+50
$\begingroup$

I have had a look at your problem. Starting afresh

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

Mathematica graphics 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]]}]

Mathematica graphics

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"}]]

Mathematica graphics

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"}]

Mathematica graphics

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.

Hopefully that helps to answer your first question.

With regard to your second question the answer depends on what you are trying to do. If you just wish to remove one frequency then a notch filter is best as I illustrated above. If you require a low pass filter then you will have to decrease your specification to avoid the numerical instability. Your values of ap and as are very difficult to implement without a high order equation. Could these be relaxed?

$\endgroup$
  • $\begingroup$ Er… so, it's impossible to design a lowpass Butterworth filter to take away the 20 Hz signal in this case? $\endgroup$ – xzczd Jun 21 '18 at 17:41
  • $\begingroup$ @xzczd This is a numerical problem so it might be possible to use longer reals. However if you are just trying to remove one frequency then a notch filter is the standard method. If you really want a low pass filter then we could search for a less restrictive specification. $\endgroup$ – Hugh Jun 21 '18 at 17:50
  • $\begingroup$ Hmm… I'm more interested in why the filter fails, I think. If that's too hard to answer, I'd like to know if there's a systematic way to decide the proper parameters for the filter model. $\endgroup$ – xzczd Jun 24 '18 at 15:02
  • $\begingroup$ I've lately been trying to simulate a Bessel filter in CUDA, using ToDiscreteTimeModel to obtain coefficients for the CUDA code from a transfer function model of an analog circuit. Unfortunately, it's creating a poorly conditioned matrix when the sampling rate is generous compared to the bandwidth, similar to the problem here. In my case, it isn't too bad: doing the calculation is single precision doesn't work, but double precision is OK, at the cost of processing time. I may someday research whether factoring the model into one and two pole filters helps. $\endgroup$ – John Doty Jun 24 '18 at 15:32
  • $\begingroup$ 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? $\endgroup$ – xzczd Jun 28 '18 at 11:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.