9
$\begingroup$

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]

Mathematica graphics

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.

$\endgroup$
  • 1
    $\begingroup$ The premise of your question is not really right -- there are many ways to design and apply linear filters other than the Butterworth filter. The tutorial here reference.wolfram.com/language/tutorial/… has several ways. $\endgroup$ – bill s Jan 11 '15 at 21:54
  • $\begingroup$ @bills As the title states I am trying to use the ButterworthFilterModel. Having selected this method I am then stuck with my questions. I could use other methods, as you suggest, but I am using the Butterworth method. $\endgroup$ – Hugh Jan 11 '15 at 22:02
  • 1
    $\begingroup$ If your goal is to "generate a filter and then apply it to some data", you are taking the long way. Modern filter design methods use techniques other than Butterworth to design the impulse response and then you can apply ListConvolve to do the filtering. $\endgroup$ – bill s Jan 11 '15 at 22:17
  • $\begingroup$ @bills Butterworth filters have a nice property of being maximally flat in the pass band. I need this property so this is why I am working with these filters. Can you better this property? $\endgroup$ – Hugh Jan 11 '15 at 22:21
  • 1
    $\begingroup$ @bills Is this discussion off-topic? I have specific questions about the Butterworth filter; you seem to be suggesting I should look elsewhere. I would like to work with what we have in Mathematica. However, if you can offer a solution for a filter that is maximally flat in the pass band and for which you can control the filter order please post. $\endgroup$ – Hugh Jan 12 '15 at 7:36
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

(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$.

$\endgroup$

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.