Building Butterworth Bandpass Filters

Question

In order to bandpass filter data it is necessary to use ButterworthFilterModel, followed by ToDiscreteTimeModel and RecurrenceFilter. I have two questions:

1. I wish to specify the order of the bandpass filter. However, this does not seem to be possible. How can I do this?

2. There is an option for ToDiscreteTimeModel called CriticalFrequency but there is no mention of prewarping which I understand to be necessary. Does CriticalFrequency introduce warping? What does it do? Must I warp?

Here is what I have worked out what to do so far. Comments on improvements are welcome. I will generate a filter and then apply it to some data. Finally by using Fourier transforms I will determine the filter I have developed. I start by defining the sample rate I will use and the high pass and low pass frequency values. I also define the attenuation I want at the cut-off frequencies.

 sr = 3500; (* Sample rate *)
 f1 = 500; (* high pass frequency (Hz) *)
 f2 = 1000; (* low pass frequency (Hz) *)
 g = -(1/2) 20 Log[10, 
  1/Sqrt[2]]; (* attenuation at cut-off frequencies *)

The reason for the 1/2 in the attenuation is that I am going to pass the data forwards and backwards through the filter and thus will get a doubling in the attenuation. I now define a warping function for the bilinear transformation which is the default option.

ClearAll[wf];
wf::usage = 
 "wf[f] gives the warped frequency in Hz for the desired frequency f \
  in Hz. Frequency must be scaled by sample rate thus values less than \
  0.5 expected.";
SetAttributes[wf, Listable];
wf[f_] := 2/(2 \[Pi]) Tan[(2 \[Pi] f)/2]

I also define a helper function to determine the order of the filter I have built

ClearAll[filterOrder];
filterOrder[bpf_] := Module[{num, den, n1, n2},
 num = bpf[[1, 1, 1, 1]];
 den = bpf[[1, 2, 1, 1]];
{n1, n2} = Flatten[Exponent[#, Variables[#]] & /@ {num, den}];
{n1, n2 - n1}
 ]

Now I wish to build a filter with high pass and low pass orders of, for example, 2. I have to supply values for the two stop frequencies and the two pass frequencies together with attention at the pass and stop frequencies. I use nominal stop frequencies and the two pass frequencies defined above. I use a value x which I vary until I get the orders I wish. How can I automatically find x?

x = 1.5;(* guess; adjust to get target filter orders *)
bpf = ButterworthFilterModel[{"Bandpass", 
 2 \[Pi] {0.9 f1, f1, f2, 1.1 f2}, {g + x, g}}];
filterOrder[bpf]

Above I have not used warping. For comparison purposes I repeat the build by first warping the target frequencies

x = 1.5;(* guess; adjust to get target filter orders *)
wbpf = ButterworthFilterModel[{"Bandpass", 
sr 2 \[Pi] wf[{0.9 f1, f1, f2, 1.1 f2}/sr], {g + x, g}}];
filterOrder[wbpf]

I now make the time domain models using ToDiscreteTimeModel. Here I could have used the option CriticalFrequency but I don't know what it does.

fm = ToDiscreteTimeModel[bpf, 1/sr];
wfm = ToDiscreteTimeModel[wbpf, 1/sr];

I would like to put all the above into one module. However because I have to adjust the attenuation x by hand this is not possible. In order to test the filters we make some data that is rich in frequencies.

nn = 100000; (* number of points *)
zv = 250; (* zero values at start and end *) 
data = 
Join[ConstantArray[0, zv], RandomReal[{-1, 1}, nn - 2 zv], 
ConstantArray[0, zv]];

The zero values are to allow for the ringing of the filter. The results of the filtering are

filterOutput = 
RecurrenceFilter[fm, Reverse[RecurrenceFilter[fm, Reverse[data]]]];
wfilterOutput = 
RecurrenceFilter[wfm, 
Reverse[RecurrenceFilter[wfm, Reverse[data]]]];

The actual filter that has been generated can be determined by calculating the Fourier transform of the output divided by the Fourier transform of the input.

fInput = Fourier[data];
fOutput = Fourier[filterOutput];
wfOutput = Fourier[wfilterOutput];
h = fOutput/fInput;
wh = wfOutput/fInput;
freqs = Table[(n - 1) sr/Length[h], {n, Length[h]}];

We can now make comparisons

p1 = Plot[{Abs[bpf[2 \[Pi] f I] bpf[2 \[Pi] f I]], 
 Abs[wbpf[2 \[Pi] f I ] wbpf[2 \[Pi] f I ]]}, {f, 0, 1500}, 
 PlotRange -> All, Frame -> True, Axes -> False, 
 FrameLabel -> {"Frequency / Hz", "Modulus"}, 
 PlotStyle -> {Blue, Brown},
 PlotLegends -> LineLegend[{"Target", "Target Warped"}],
 Epilog -> {Pink, Line[{{f1, 0}, {f1, 1.1}}], 
 Line[{{f2, 0}, {f2, 1.1}}]}];
 p2 = ListPlot[{Transpose[{freqs, Abs[h]}][[3000 ;; nn/2]], 
 Transpose[{freqs, Abs[wh]}][[3000 ;; nn/2]]}, Joined -> True, 
  Frame -> True, PlotStyle -> {Green, Orange},
  PlotLegends -> LineLegend[{"Calculated", "Calculated Warped"}]];
 Show[p1, p2]

The pink lines are at the cut-off frequencies. It can be seen that the spectrum calculated using the warped filter (orange) agrees best with the target filter (blue). Without the warping (green) the calculated filter is badly off. The warped target spectrum (brown) looks wrong but is the one to use.

A possibility for calculating the filter is to use a version of FindRoot. However, the results are integers so this does not work.

Thanks for help in answering the questions.

Answer 1

As of version 10 (certainly in version 10.3) you can specify the filter order and the band pass and band stop frequencies. This is what I wanted in my first question.

I cannot find any reference to prewarping. I guess you just have to know that this is necessary.

Answer 2

(to answer the second part)

"CriticalFrequency" is a suboption to the "BilinearTransform" (aka Tustin's) method. Bilinear transformation introduces warping, because it maps the entire imaginary axis onto the circumference of the unit circle. With "CriticalFrequency"->$f$, there is no distortion at frequency $f$.