2
$\begingroup$

I work in Mathematica with a Butterworth filter, the transfer function of which depends on the selected order and cutoff frequency.

I want to evaluate the change in the location of the roots on the complex plane if I slightly modify the frequency response of the filter. To do this, I create a data set of the amplitude-frequency characteristics, select a model with custom parameters (the roots are custom parameters) and fit it to the data using FindFit.

There are no problems with this, but the problem in another thing.

If I choose the high order of the model ($n=8$ and above), then the speed of the FindFit will decrease sharply. For example, for the order $n = 8$ it took me about 15 - 20 minutes.

I tried to change the calculation method, but it didn’t really affect anything.

I try different methods, and the usual FindFit, and Nminimize.

Are there any tricks to increase the calculation speed and make it more or less convenient?

Here is my code:

ClearAll["Global`*"]

pars = {n = 3} (*Order of model*)

btf = N[First@ButterworthFilterModel[{"Lowpass", 3, 1}, s]]

tf = btf[[1, 1, 1]]/btf[[2, 1, 1]];

s = ω I;

r = ComplexExpand[Re[tf]] // Simplify;

y = ComplexExpand[Im[tf]] // Simplify;

afc[ω_] = 
 Sqrt[r^2 + y^2] // Simplify;(*Amplitude-Frequency Characteristic*)

data = {#, afc[#]} & /@ Table[ω, {ω, 0, 2, .05}];

lp = ListPlot[data, PlotStyle -> {PointSize[0.01]}, 
   DisplayFunction -> Identity];

Show[lp, DisplayFunction -> $DisplayFunction, 
 PlotRange -> Full]

pol = 1/Product[(s - (Subscript[a, k] + Subscript[b, k] I)), {k, 1, n}]

R = ComplexExpand[Re[pol]] // Simplify;

Y = ComplexExpand[Im[pol]] // Simplify;

parameters = Join[Table[Subscript[a, k], {k, 1, n}], Table[Subscript[b, k], {k, 1, n}]]

constraints = Table[Sign[Subscript[a, k]] <= 0, {k, 1, n}]

model = Sqrt[R^2 + Y^2] // Simplify

f1 = FindFit[data, model, parameters, ω]

f2 = FindFit[data, {model, constraints}, parameters, ω, Method -> NMinimize]

p1 = Show[lp, DisplayFunction -> $DisplayFunction, PlotRange -> Full];

p2 = LogLinearPlot[Evaluate[model /. {f1, f2}], {ω, 0, 10}, 
   PlotStyle -> {Thick, Dashed}, PlotRange -> Full];

GraphicsRow[{p1, p2}]

EDIT:

I replaced the lines in the code:

R = ComplexExpand [Re [pol]] // Simplify; 

Y = ComplexExpand [Im [pol]] // Simplify; 

model = Sqrt [R ^ 2 + Y ^ 2] // Simplify 

On the:

model = Abs [pol] // Simplify; 

It began to work faster.

$\endgroup$
9
  • 1
    $\begingroup$ The part that hangs is not the FindFit, but the ComplexExpand and Simplify for both R and Y. You can replace these with Re[pol], Im[pol], don't do the expansion or simplify and it completes and I get a plot. $\endgroup$ – flinty Jun 17 '20 at 14:00
  • $\begingroup$ I replaced the lines in the code R = ComplexExpand [Re [pol]] // Simplify; Y = ComplexExpand [Im [pol]] // Simplify; model = Sqrt [R ^ 2 + Y ^ 2] // Simplify On the: model = Abs [pol] // Simplify; It began to work faster. $\endgroup$ – dtn Jun 17 '20 at 14:08
  • $\begingroup$ And did it work? $\endgroup$ – flinty Jun 17 '20 at 14:10
  • $\begingroup$ It began to work faster. $\endgroup$ – dtn Jun 17 '20 at 14:10
  • 1
    $\begingroup$ I noticed you're not using n when you set up your filter. Also I think you should have a look at TransferFunctionPoles and TransferFunctionZeros. e.g TransferFunctionPoles[ButterworthFilterModel[{"Lowpass", n, 1}, s]] - and it doesn't have any zeros, only poles. $\endgroup$ – flinty Jun 17 '20 at 16:25
2
$\begingroup$

You can get the poles of the Butterworth filter with different lowpass filter cutoffs without FindFit. Let me know if this is helpful:

order = 24;
cutoff1 = 1.5;
cutoff2 = 2.4;
btf1 = ButterworthFilterModel[{"Lowpass", order, cutoff1}, s];
btf2 = ButterworthFilterModel[{"Lowpass", order, cutoff2}, s];
poles1 = Flatten@TransferFunctionPoles[btf1];
poles2 = Flatten@TransferFunctionPoles[btf2];
topoints[complexes_] := Through[{Re, Im}[#]] & /@ complexes;
Plot[{Abs[btf1[I f]], Abs[btf2[I f]]}, {f, 0, 4}, 
 AxesLabel -> {"f", "α"}, 
 PlotLegends -> {"c=" <> ToString@cutoff1, "c=" <> ToString@cutoff2}]
Graphics[{
  Red,    Point[topoints[poles1]],
  Orange, Point[topoints[poles2]]
}, Axes -> True, Frame -> True, PlotRange -> All]

enter image description here

$\endgroup$
6
  • $\begingroup$ This can be used to get the roots and use them as initial parameter values. So, the answer is useful, but it is better to pay attention to the reasons for a sharp drop in the speed of calculations when working with a large number of optimized parameters. $\endgroup$ – dtn Jun 17 '20 at 19:16
  • $\begingroup$ Could you elaborate on what you're actually doing to the roots, mathematically? I got the impression you wanted to study the change in the roots when you change the cutoff, which is what this shows - in which case FindFit is not needed as it's too costly. $\endgroup$ – flinty Jun 17 '20 at 19:22
  • $\begingroup$ Yes, I’ll try to clarify. You change the parameters of the already known Butterworth filter structure. And I need to evaluate the location of the roots on the complex plane, if I change its frequency response, for example, add a slope of -20 dB / dec to the cutoff frequency. $\endgroup$ – dtn Jun 17 '20 at 19:30
  • 1
    $\begingroup$ If something else needs to be clarified, then I will be happy to answer you. $\endgroup$ – dtn Jun 17 '20 at 19:30
  • 2
    $\begingroup$ Presumably that's this part of your code: tf = btf[[1, 1, 1]]/btf[[2, 1, 1]]. If you turn this back into a TransferFunctionModel you can use my approach above without needing to use FindFit. There is also TransferFunctionTransform if you want to tweak an existing TransferFunctionModel - and this is recommended vs. modifying the internals. $\endgroup$ – flinty Jun 17 '20 at 19:37

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.