diff --git a/scripts/plot.sh b/scripts/plot.sh index 4cf8767..1c16c39 100755 --- a/scripts/plot.sh +++ b/scripts/plot.sh @@ -27,11 +27,12 @@ case $1 in "lags"|"lag"|"delay"|"delays") gnuplot_file=timelag_atlas.gp gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g') - scripts/propagate_tables.sh - for tabfile in data/tables/timelag_*${errtype}*.tab; + for tabfile in data/tables/lag_*.tab; do - ref_band_extracted=$(basename $tabfile|sed 's|timelag_\([^≺]*\)[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*}.tab|\1|') - echo_band=$(basename $tabfile|sed 's|timelag_[^≺]*[_ ]≺[_ ]\([^≺_ ]*\)[_ ]{[^_ ]*}.tab|\1|') + ref_band_extracted=$(basename $tabfile| + sed 's|lag_\([0-9]\{4\}A\)_[0-9]\{4\}A.tab|\1|') + echo_band=$(basename $tabfile| + sed 's|lag_[0-9]\{4\}A_\([0-9]\{4\}A\).tab|\1|') if [[ "$echo_band" == "$ref_band" ]] ; then continue; fi gnuplot_input_edit=$(echo "$gnuplot_input"| sed "s|%FILE%|$tabfile|"| diff --git a/scripts/tophat_fft.pl b/scripts/tophat_fft.pl new file mode 100755 index 0000000..c2f2979 --- /dev/null +++ b/scripts/tophat_fft.pl @@ -0,0 +1,85 @@ +#!/usr/bin/env perl + +use feature say; +use utf8; +use PDL; +use PDL::FFT; +use PDL::IO::Misc; + +my $_2π = 2*3.1415926539; + +# space parameters +our $xres= 0.1; +our $xmin = 0; +our $xmax = 2500; +our $xbins = ($xmax-$xmin)/$xres + 1; + +say "Using $xbins bins"; + +# tophat parameters +my $μ1=7; +my $Δ1=5; # full width + +my @tophat_list = ([$μ1,$Δ1],[28,8],[15,11]); + +our $tophat_count = 0; +foreach (@tophat_list) { + $tophat_count++; + + # Capture tophat into piddle + # Complex Number z(x) = u(x) + iv(x) + my $u = tophat(@$_[0],@$_[1]); + my $v = zeroes($u); + + $x_coords = $u->xlinvals($xmin,$xmax); + + # Output x-domain tophat + wcols $x_coords,$u,"analyses/tables/tophat${tophat_count}.tab"; + + # Transform z(x) to Z(x⁻¹) = U(x⁻¹) + iV(x⁻¹) + my $U = $u; + my $V = $v; + fft($U,$V); + + # Determine frequency coordinates + my $num_elements = nelem($U); + say "Found $num_elements elements."; + + $f = ($num_elements % 2 == 0) ? + $U->xlinvals( + -(${num_elements}/2-1)/${num_elements}/${xres}, + 1/2/${xres} + )->rotate(-(${num_elements}/2 -1)) + : + $U->xlinvals( + -(${num_elements}/2-0.5)/${num_elements}/${xres}, + (${num_elements}/2-0.5)/${num_elements}/${xres} + )->rotate(-(${num_elements}-1)/2); + + # Output frequency-domain tophat + wcols $f, $U, $V, "fft${tophat_count}.tab"; + + # Calculate x offset + # This currently multiplies the imaginary component by -1, + # and I really need to figure out why this is necessary for + # proper output. + my $φdiff = atan2(-1*$V,$U); + + wcols $f,$φdiff,"analyses/tables/tophat_φdiff${tophat_count}.tab"; + + my $offset = $φdiff/($_2π*$f); + + # Output frequency-domain time delay for given tophat + wcols $f,$offset,"analyses/tables/tophat_fft${tophat_count}.tab"; +} + +sub tophat { + (my $mean,my $width) = @_; + my $halfwidth = $width/2; + my @vals = (); + for (my $x_coord = $xmin; $x_coord <= $xmax; $x_coord += $xres) { + push @vals, ($x_coord >= ($mean - $halfwidth) && + $x_coord <= ($mean + $halfwidth)) ? 1/$width : 0; + } + return pdl(@vals); +} \ No newline at end of file