diff --git a/plot.sh b/plot.sh index d4ceb99..9b4b55c 100755 --- a/plot.sh +++ b/plot.sh @@ -30,7 +30,6 @@ case $1 in done echo "$gnuplot_input"|perl -pe 's|␤|\n|g' > ${gnuplot_file} gnuplot $gnuplot_file - ;; "lags"|"lag"|"delay"|"delays") @@ -52,10 +51,10 @@ case $1 in ;; "tophat") - gnuplot_file=tophat_w_fft.gp + gnuplot_file="tophat_w_fft.gp" ;; - + *) echo "Did not understand plot type." - ;; + ;; esac \ No newline at end of file diff --git a/scripts/tophat_fft.pl b/scripts/tophat_fft.pl index de9faa0..a1821ea 100755 --- a/scripts/tophat_fft.pl +++ b/scripts/tophat_fft.pl @@ -1,7 +1,7 @@ #!/usr/bin/env perl use feature say; -use utf8; +use utf8; use PDL; use PDL::FFT; use PDL::IO::Misc; @@ -9,38 +9,52 @@ use PDL::IO::Misc; my $twoPI = 2*3.1415926539; # space parameters -our $Δt=.001; +our $xres=.001; our $xmin = 0; our $xmax = 250; -our $num_bins = ($xmax-$xmin)/$Δt; +our $num_bins = ($xmax-$xmin)/$xres; + +# tophat parameters +my $μ1=15; +my $Δ1=5; # full width +my $μ2=15; +my $Δ2=5; # full width +my $μ3=15; +my $Δ3=5; # full width say "Creating $num_bins bins"; # Complex Number z(x) = u(x) + iv(x) -my @th_list=tophat(10,7); +my @th_list=tophat($μ1,$Δ1); my $u = pdl(@th_list); my $v = zeroes($u); -# Output time-domain tophat -open $tophattab, '>', 'analyses/tables/tophat.tab'; -for ($i=0; $i < $num_bins; $i++) { - my $x_coord = ($xmax-$xmin)*($i/$num_bins) + $xmin; - say $tophattab "$x_coord $th_list[$i]"; + +foreach (('1','2','3')) { + + # Output time-domain tophat + open $tophattab, '>', "analyses/tables/tophat${$_}.tab"; + for ($i=0; $i < $num_bins; $i++) { + my $x_coord = ($xmax-$xmin)*($i/$num_bins) + $xmin; + say $tophattab "$x_coord $th_list[$i]"; + } + close $tophattab; + + # Transform z(x) to Z(1/x) = U(1/x) + iV(1/x) + my $U = $u; + my $V = $v; + fft($U,$V); + # Determine frequency coordinates + my $num_elements = nelem($U); + say "Found $num_elements elements."; + $f = $U->xlinvals(-($num_elements/2-1)/$num_elements/$xres,1/2/$xres)->rotate(-($num_elements/2 -1)); + + my $φdiff = atan2($V,$U); + my $timelag = $φdiff/($twoPI*$f); + + wcols $f,$timelag,"analyses/tables/tophat_fft${$_}.tab"; + } -close $tophattab; - -# Transform z(x) to z(1/x) = u(1/x) + iv(1/x) -fft($u,$v); -# Determine frequency coordinates -my $num_elements = nelem($u); -say "Found $num_elements elements."; -$f = $u->xlinvals(-($num_elements/2-1)/$num_elements/$Δt,1/2/$Δt)->rotate(-($num_elements/2 -1)); - -my $φdiff = atan2($v,$u); -my $timelag = $φdiff/($twoPI*$f); - -wcols $f,$timelag,'analyses/tables/tophat_fft.tab'; - sub tophat { (my $mean,my $width) = @_;