Menu
▾
▴
Re: Adding unwrap to gnuplot
|
From: Jonathan T. <jt...@as...> - 2014-01-30 18:12:02
|
On Thu, Jan 30, 2014 at 02:00:18PM +0000, Jon Gjengset wrote:
> When plotting phase data in gnuplot, you often encounter "jumps" in a
> graph when a value wraps around 2??. A common way of dealing with this is
> to apply an unwrap[1] function. What this function does is add or remove
> 2pi from the next value to be plotted (v) as long as the difference
> between v and the previous value (w) is outside the range [-pi,pi].
[[...]]
Here is a perl script which does such an "unwrapping"; feel free to
use it. (I realise that this isn't quite what you asked for, which
was a gnuplot builtin function. But my perl script might be a useful
substitute for some use cases.)
--- begin perl script ---
#!/usr/bin/perl -w
# $Header: /home/jonathan/CVSROOT/src/misc/unwrap.phase,v 1.3 2006/12/19 14:33:47 jonathan Exp $
################################################################################
#
# unwrap.phase -- unwrap the phase of a time series of angles
# This program is copyright (C) 2006 by Jonathan Thornburg <jt...@ae...>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
#
################################################################################
my $pi = 4.0*atan2(1.0,1.0);
my $default_period = 2.0*$pi;
my $default_input_column = 1;
################################################################################
my $help_msg = <<"EOF"; # expand variables
This program "unwraps" the phase of a time series of "angles". That is,
given a time series, this function adds an integer multiple of a specified
period to each value so as to minimize the total of the absolute differences
of the sums,
output_i := input_i + N_i*period
N_i chosen to minimize |output_i - output_{i-1}|
Usage:
unwrap.phase [ --period <number> ]
[ --input-column <integer> ]
[ --debug <integer> ]
<input >output
or
unwrap.phase --help # print this output
or
unwrap.phase --help-example # print this output
Any all-whitespace or comment lines in the input ('#' in column 1) are
passed through to the output unchanged. Otherwise, each input line is
split into whitespace-separated fields, and an angle is extracted. The
corresponding output line is a copy of the input line, followed by N_i
and output_i (tab-separated).
The --period optional argument (default 2*pi) specifies the period.
The --input-column optional argument (default ${default_input_column}) specifies the 1-origin
column containing the input angle.
EOF
################################################################################
my $help_example_msg = <<'EOF';
% cat unwrap.phase.test-data
# column 1 = input
# column 2 = desired output
0.0 0.0
0.3 0.3
0.6 0.6
0.9 0.9
0.2 1.2
0.5 1.5
0.8 1.8
0.1 2.1
0.4 2.4
0.7 2.7
0.0 3.0
0.9 2.9
0.5 2.5
0.1 2.1
0.7 1.7
0.3 1.3
0.9 0.9
0.5 0.5
0.1 0.1
0.7 -0.3
0.3 -0.7
0.9 -1.1
0.1 -0.9
% unwrap.phase --period 1 <unwrap.phase.test-data
# column 1 = input
# column 2 = desired output
0.0 0.0 0 0
0.3 0.3 0 0.3
0.6 0.6 0 0.6
0.9 0.9 0 0.9
0.2 1.2 1 1.2
0.5 1.5 1 1.5
0.8 1.8 1 1.8
0.1 2.1 2 2.1
0.4 2.4 2 2.4
0.7 2.7 2 2.7
0.0 3.0 3 3
0.9 2.9 2 2.9
0.5 2.5 2 2.5
0.1 2.1 2 2.1
0.7 1.7 1 1.7
0.3 1.3 1 1.3
0.9 0.9 0 0.9
0.5 0.5 0 0.5
0.1 0.1 0 0.1
0.7 -0.3 -1 -0.3
0.3 -0.7 -1 -0.7
0.9 -1.1 -2 -1.1
0.1 -0.9 -1 -0.9
EOF
################################################################################
use strict;
use Getopt::Long;
my $true = 1;
my $false = 0;
my $fuzz = 1.0e-10;
my $period = $default_period;
my $input_column = $default_input_column;
my $debug = 0;
my $help_flag = $false;
my $help_example_flag = $false;
$Getopt::Long::autoabbrev = $false; # forbid abbreviations
$Getopt::Long::getopt_compat = $false; # forbid starting options with "+"
GetOptions(
'period=f' => \$period,
'input-column=i' => \$input_column,
'debug=i' => \$debug,
'help' => \$help_flag,
'help-example' => \$help_example_flag,
)
|| die $help_msg;
if ($help_flag)
{ print $help_msg; exit 0; }
if ($help_example_flag)
{ print $help_example_msg; exit 0; }
if ($period == 0.0)
{ die "period must be nonzero!\n"; }
if ($debug > 0)
{
print "using period=${period}\n";
}
if ($debug == 99)
{
for (my $x = -1.6 ; fuzzy_LE($x,1.6) ; $x += 0.1)
{
print $x, "\t", round_to_nearest_int($x), "\n";
}
exit 0;
}
# perl uses 0-origin fields internally
$input_column -= 1;
my $previous_output = undef();
my $input_line_number = 0;
while (my $line = <STDIN>)
{
++$input_line_number;
if (($line =~ /^\s$/) || ($line =~ /^#/))
{ print $line; next; } # *** LOOP CONTROL ***
chomp $line;
my @fields = split(/\s+/, $line);
my $input = $fields[$input_column];
if (! defined($input))
{ die "couldn't find input field in input line ${input_line_number}!\n"; }
$input += 0.0; # make sure it's numeric
my $N;
if (defined($previous_output))
{
# choose $N to minimize absolute difference
# |$output - $previous_output|
# where $output is computed below
my $Nf = ($previous_output - $input) / $period;
$N = round_to_nearest_int($Nf);
if ($debug > 5)
{ print "input=${input} ==> Nf=${Nf} N=${N}\n"; }
}
else {
$N = 0;
}
if ($debug > 0)
{ print "input=${input} ==> N=${N}\n"; }
my $output = $input + $N*$period;
print $line, "\t", $N, "\t", $output, "\n";
$previous_output = $output;
}
################################################################################
#
# This function rounds a floating-point number to the nearest integer.
#
# Arguments:
# $x = The number to round.
#
sub round_to_nearest_int
{
my ($x) = @_;
# remaining code only handles values >= 0
if ($x < 0) { return -round_to_nearest_int(-$x); }
return int(0.5+$x);
}
################################################################################
#
# These functions do fuzzy comparisons of real numbers, returning
# Boolean values $true or $false.
#
# Global variables:
# $fuzz
#
sub fuzzy_EQ
{
my ($x,$y) = @_;
return (abs($x-$y) <= $fuzz) ? $true : $false;
}
sub fuzzy_LE
{
my ($x,$y) = @_;
return (($x < $y) || fuzzy_EQ($x,$y)) ? $true : $false;
}
--- end perl script ---
ciao,
--
-- "Jonathan Thornburg [remove -animal to reply]" <jt...@as...>
Dept of Astronomy & IUCSS, Indiana University, Bloomington, Indiana, USA
"There was of course no way of knowing whether you were being watched
at any given moment. How often, or on what system, the Thought Police
plugged in on any individual wire was guesswork. It was even conceivable
that they watched everybody all the time." -- George Orwell, "1984"
|
