Science/psdlag-agn
1
0
Fork 0
mirror of https://asciireactor.com/otho/psdlag-agn.git synced 2024-11-24 12:35:06 +00:00
Code Issues Packages Projects Releases Wiki Activity

tophat transfer function basically done but for some reason proper output requires V(x⁻¹) to be multiplied by ι²

Browse Source
This commit is contained in:
caes 2016-08-14 01:41:32 -04:00
parent 6ffe2c92bc
commit 319a5766ca
4 changed files with 25 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
plot.sh
View File

@ -10,13 +10,10 @@ echo_bands=$(ls analyses/*CM*|sed 's|[^≺]*≺_\(.\{5\}\).*|\1|')
# echo Using list of reverberated bands: $echo_bands # echo Using list of reverberated bands: $echo_bands
echo Propagating tables.
scripts/propagate_tables.sh > /dev/null
case $1 in case $1 in
"PSD"|"psd"|"PSDs"|"PSDS"|"psds") "PSD"|"psd"|"PSDs"|"PSDS"|"psds")
echo "Producing PSD atlas."
gnuplot_file=psd_atlas.gp gnuplot_file=psd_atlas.gp
scripts/propagate_tables.sh
gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g') gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g')
for tabfile in analyses/tables/PSD_*${errtype}*.tab; for tabfile in analyses/tables/PSD_*${errtype}*.tab;
do do
@ -33,9 +30,9 @@ case $1 in
;; ;;
"lags"|"lag"|"delay"|"delays") "lags"|"lag"|"delay"|"delays")
echo "Producing time delay atlas."
gnuplot_file=timelag_atlas.gp gnuplot_file=timelag_atlas.gp
gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g') gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g')
scripts/propagate_tables.sh
for tabfile in analyses/tables/timelag_*${errtype}*.tab; for tabfile in analyses/tables/timelag_*${errtype}*.tab;
do do
ref_band_extracted=$(basename $tabfile|sed 's|timelag_\([^≺]*\)[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*}.tab|\1|') ref_band_extracted=$(basename $tabfile|sed 's|timelag_\([^≺]*\)[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*}.tab|\1|')
@ -51,7 +48,9 @@ case $1 in
;; ;;
"tophat") "tophat")
gnuplot_file="tophat_w_fft.gp" scripts/tophat_fft.pl
gnuplot scripts/templates/tophat_freqdomain.gp
gnuplot scripts/templates/tophat_timedomain.gp
;; ;;
*) *)

BIN
pres/img/tophat_freqdomain.pdf
View File

Binary file not shown.

BIN
pres/img/tophat_timedomain.pdf
View File

Binary file not shown.

38
scripts/tophat_fft.pl
View File

@ -23,13 +23,12 @@ my $Δ1=5; # full width
my @tophat_list = ([$μ1,$Δ1],[28,8],[15,11]); my @tophat_list = ([$μ1,$Δ1],[28,8],[15,11]);
our $tophat_count = 0; our $tophat_count = 0;
foreach my $pars (@tophat_list) { foreach (@tophat_list) {
$tophat_count++; $tophat_count++;
# Capture tophat into piddle # Capture tophat into piddle
# Complex Number z(x) = u(x) + iv(x) # Complex Number z(x) = u(x) + iv(x)
my @tophat=tophat(@$pars[0],@$pars[1]); my $u = tophat(@$_[0],@$_[1]);
my $u = pdl(@tophat);
my $v = zeroes($u); my $v = zeroes($u);
$x_coords = $u->xlinvals($xmin,$xmax); $x_coords = $u->xlinvals($xmin,$xmax);
@ -46,25 +45,28 @@ foreach my $pars (@tophat_list) {
my $num_elements = nelem($U); my $num_elements = nelem($U);
say "Found $num_elements elements."; say "Found $num_elements elements.";
if ($num_elements % 2 == 0) { $f = ($num_elements % 2 == 0) ?
# even number of bins $U->xlinvals(
$f = $U->xlinvals( -(${num_elements}/2-1)/${num_elements}/${xres},
-(${num_elements}/2-1)/${num_elements}/${xres}, 1/2/${xres} 1/2/${xres}
)->rotate(-(${num_elements}/2 -1)); )->rotate(-(${num_elements}/2 -1))
} :
else { $U->xlinvals(
#odd number of bins -(${num_elements}/2-0.5)/${num_elements}/${xres},
$f = $U->xlinvals( (${num_elements}/2-0.5)/${num_elements}/${xres}
-(${num_elements}/2-0.5)/${num_elements}/${xres}, )->rotate(-(${num_elements}-1)/2);
(${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 # Calculate x offset
my $φdiff = atan2($V,$U); # 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(-$V,$U);
my $offset = $φdiff/($_2π*$f); my $offset = $φdiff/($_2π*$f);
# Output frequency-domain time delay for given tophat
wcols $f,$offset,"analyses/tables/tophat_fft${tophat_count}.tab"; wcols $f,$offset,"analyses/tables/tophat_fft${tophat_count}.tab";
} }
@ -77,5 +79,5 @@ sub tophat {
push @vals, ($x_coord >= ($mean - $halfwidth) && push @vals, ($x_coord >= ($mean - $halfwidth) &&
$x_coord <= ($mean + $halfwidth)) ? 1/$width : 0; $x_coord <= ($mean + $halfwidth)) ? 1/$width : 0;
} }
return @vals; return pdl(@vals);
} }