bode plots
Michiel Kreutzer
I am trying to make a couple of Bode plots using gnuplot.
For that I need to plot the atan(Im(f)/Re(f)), where
f is a complex function.
For several functions, f = exp(-iw) to name just one,
the function starts to oscilate very fast from pi to -pi.
Of course, the atan of a number is strictly speaking
not between -pi and pi, but one can add or substract
2*k*pi from it, where k is any integer.
If I could find a way to increment k by one each time
the function jumps from -pi to pi, I would get a
smooth line.
For that I need to have the derivative, which I could get
from :
f(x) = atan(...)
eps = 0.00001
dfdx = (f(x+eps) - f(x)) / eps
If the derivative is positive (indicating the jump),
k has to be incremented by one, and then
g(x) = f(x) - 2*k*pi
would give the smooth plot.
does anybody know a way to do this in gnuplot?
I have added a small gnuplot file to illustrate my problem.
cheers
Michiel
----------------------------
f(x) = exp({0,1}*x*-2*pi)
set terminal pslatex rotate auxfile
set output "bode_pfr.tex"
set size 1,1.6
set origin 0.1, 0.1
set multiplot
set size 0.8,0.8
set origin 0.1,0.9
set lmargin 12
set border 15 lw 0.5
set nokey
set logscale
set format y "$10^{%T}$"
set format x "$10^{%T}$"
set title
set ylabel "$|G(i \\omega)|$"
set xlabel "$f$"
plot [x=0.01:10] [1e-4:1e2] abs(f(x))
set size 0.8,0.8
set origin 0.1,0.1
set lmargin 12
set nologscale y
set format y "$%2.1P \\pi$"
set ylabel "arg $(G(i \\omega))$"
set key left
set ytic -1.5*3.141593,0.5*3.14193,1*3.141593
plot [x=0.01:10] [-5:3.2] arg(f(x))
set nomultiplot
set output
------------------------------------------------------------
Hans-Bernhard Broeker
> For that I need to plot the atan(Im(f)/Re(f)), where
> f is a complex function.
It's a lot easier to write arg(f) instead. Or, if you really want to
use atan, use atan2, instead:
atan2(imag(f), real(f))
Atan is wrong because it only has a range of -pi/2 to pi/2, i.e. it
will return the same angle for atan(x/y) and atan(-x/-y). Atan2 gets
that right.
> For several functions, f = exp(-iw) to name just one,
> the function starts to oscilate very fast from pi to -pi.
Once you've called arg() or atan2() on the result, there's no reliable
way you can fix what this broke. If you want the "real" argument,
whatever that may mean, you'll have to express the function itself in
terms of radius and argument. In the case at hand, that's simple:
real_angle_of_f(x) = -2*pi*x
> If I could find a way to increment k by one each time
> the function jumps from -pi to pi, I would get a
> smooth line.
In the general case, that's impossible. You'll never quite be able to
tell the difference between an actual jump (i.e.a discontinuity) and a
place where the function just changes rapidly, just by looking at a
finite set of computed function values.
> does anybody know a way to do this in gnuplot?
Don't. Plot in polar coordinates instead, and you'll be a lot better
off.
set polar
set parametric
plot arg(f(t)), t
The real hard part of the original problem would be that you would
have to check not only the position currently being plotted, but
rather the *whole*x*range* up to there to check for any possible jumps
that have to be compensated.
--
Hans-Bernhard Broeker (bro...@physik.rwth-aachen.de)
Even if all the snow were burnt, ashes would remain.
Michiel Kreutzer
> > For that I need to plot the atan(Im(f)/Re(f)), where
> > f is a complex function.
>
> It's a lot easier to write arg(f) instead. Or, if you really want to
> use atan, use atan2, instead:
I fully agree, and -- if you look at the example I have provided
in the original post -- you will see that that is what I do, ie:
plot [x=0.01:10] [-5:3.2] arg(f(x))
> Once you've called arg() or atan2() on the result, there's no reliable
> way you can fix what this broke. If you want the "real" argument,
Not "real", but rather more informative for a review of my topic
(see below).
> whatever that may mean, you'll have to express the function itself in
> terms of radius and argument. In the case at hand, that's simple:
>
> real_angle_of_f(x) = -2*pi*x
Yep, but this is a trivial one, and for the G(iw) = exp(-iw) case,
that is what I did to get the "right" plot. If I could have solved
my problems with elementary calculus, I would have ;-). But here is
a different one for you:
q(x,Bo) = sqrt(1+{0,1}*2*pi*x*4/Bo)
f(x,Bo) = 4*q(x,Bo)/ \
((1+q(x,Bo))**2*exp(Bo/(-2)*(1-q(x,Bo))) - \
(1-q(x,Bo))**2*exp(Bo/(-2)*(1+q(x,Bo))))
Now the simple way out does not apply anymore! But what bothers
me is that the Bode plots, in for instance the textbooks where
I picked up the aforementioned one, they all have these very
nice "continuous" bode plots.
> In the general case, that's impossible. You'll never quite be able to
> tell the difference between an actual jump (i.e.a discontinuity) and a
> place where the function just changes rapidly, just by looking at a
> finite set of computed function values.
>
> > does anybody know a way to do this in gnuplot?
>
> Don't. Plot in polar coordinates instead, and you'll be a lot better
> off.
I'll forward your suggestion to my thesis supervisor. However, in what
I am doing now (deconvolution of retention time distribution using
FFT), I think the Bode plots are quite often a "easier" fingerprint
of the situation at hand than the polar plots (I think those are
sometines referred to as Nyquist plots, but my math-textbook is
at the lab now, so I am not sure)
Anyway, I am getting off-topic now: I think the
best way out is just writing a do-loop in a simple FORTRAN program
and incrementing k at each jump. That will definitely work (also
because I know that the argument of my functions allways have a
negative derivative - Q not ED ). I was just hoping it could be
done inside gnuplot, because I have to do quite a lot of these!
Thanks anyway for your very quick reply!
cheers
Michiel
Michiel Kreutzer
> > For that I need to plot the atan(Im(f)/Re(f)), where
> > f is a complex function.
>
> It's a lot easier to write arg(f) instead. Or, if you really want to
> use atan, use atan2, instead:
I fully agree, and -- if you look at the example I have provided
in the original post -- you will see that that is what I do, ie:
plot [x=0.01:10] [-5:3.2] arg(f(x))
> Once you've called arg() or atan2() on the result, there's no reliable
> way you can fix what this broke. If you want the "real" argument,
Not "real", but rather more informative for a review of my topic
(see below).
> whatever that may mean, you'll have to express the function itself in
> terms of radius and argument. In the case at hand, that's simple:
>
> real_angle_of_f(x) = -2*pi*x
Yep, but this is a trivial one, and for the G(iw) = exp(-iw) case,
that is what I did to get the "right" plot. If I could have solved
my problems with elementary calculus, I would have ;-). But here is
a different one for you:
q(x,Bo) = sqrt(1+{0,1}*2*pi*x*4/Bo)
f(x,Bo) = 4*q(x,Bo)/ \
((1+q(x,Bo))**2*exp(Bo/(-2)*(1-q(x,Bo))) - \
(1-q(x,Bo))**2*exp(Bo/(-2)*(1+q(x,Bo))))
Now the simple way out does not apply anymore! But what bothers
me is that the Bode plots, in for instance the textbooks where
I picked up the aforementioned one, they all have these very
nice "continuous" bode plots.
> In the general case, that's impossible. You'll never quite be able to
> tell the difference between an actual jump (i.e.a discontinuity) and a
> place where the function just changes rapidly, just by looking at a
> finite set of computed function values.
>
> > does anybody know a way to do this in gnuplot?
>
> Don't. Plot in polar coordinates instead, and you'll be a lot better
> off.
I'll forward your suggestion to my thesis supervisor. However, in what
Hans-Bernhard Broeker
> Now the simple way out does not apply anymore! But what bothers
> me is that the Bode plots, in for instance the textbooks where
> I picked up the aforementioned one, they all have these very
> nice "continuous" bode plots.
That they may do --- for the few plots needed in an expensively
produced textbook, you can always hire someone to fiddle with the step
function to "correct" the discontinuities until you like it. It's a
whole different story if you try to achieve the same result fully
automatically.
> Anyway, I am getting off-topic now: I think the
> best way out is just writing a do-loop in a simple FORTRAN program
> and incrementing k at each jump.
If you really think you have to, you _can_ do that in gnuplot, too.
But it'll be quite slow in the plotting. Use recursion and the '?'
operator. For the example of an argument known to be strictly
monotonously decreasing:
# set up some parameters
xmin = 0.0; xmax = 20.0; searchstep = 0.2
# make sure the sampling done by gnuplot matches the 'searchstep':
set xrange [ xmin : xmax] noreverse nowriteback
set samples (xmax - xmin) / searchstep + 1
# define the function
f(x) = arg(exp(-{0,1}*2*pi*x))
# and the "correction term"
stepfunc(x,x0) = (x < x0) ? 0 : \
(stepfunc(x - searchstep, x0) + \
((f(x) > f(x - searchstep)) ? -2*pi : 0))
# show that the plot actually works:
plot f(x), stepfunc(x), f(x)+stepfunc(x)
> done inside gnuplot, because I have to do quite a lot of these!
Exactly if you have to do quite a lot of these, you'll not want to do
them inside gnuplot --- it's just too damned slow.
Michiel Kreutzer
[snipped discussion "beautifying" bode plots]
> > I think the best way out is just writing a do-loop
> > in a simple FORTRAN program and incrementing k at
> > each jump.
>
> If you really think you have to, you _can_ do that in gnuplot, too.
> But it'll be quite slow in the plotting. Use recursion and the '?'
> operator. For the example of an argument known to be strictly
> monotonously decreasing:
[snipped perfect demo of the method]
Perfect, works like a charm!
> > done inside gnuplot, because I have to do quite a lot of these!
>
> Exactly if you have to do quite a lot of these, you'll not want to do
> them inside gnuplot --- it's just too damned slow.
Takes a few seconds to make the plots here. Fast enough for me. Faster
then going the FORTRAN route, anyway...
Thanks a lot for your inputs, you have been very helpful!
cheers
Michiel