#!/usr/bin/perl
#perl smoothconv.pl static.txt2 std wind col pad
#example perl smoothconv.pl static.txt 70 200 2 3
#70: standard deviation
#200: gaussian window
#2: column you would like to smooth from input
#3: padding (scale) >>> to remove edge effect
$file=$ARGV[0];
$sd=$ARGV[1];
$gw=$ARGV[2];
$col=$ARGV[3];
$pad=$ARGV[4];
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
#push @ax1, @row; #only one col input
push @ax1, \@row;
}
close( FILE );
open OUT1," > in.txt" or die "$!\n";
for my $i (0..$#ax1-1)
{
push @ax, $ax1[$i][$col-1];
push @ax_abs, abs($ax1[$i][$col-1]);
print OUT1 $ax1[$i][$col-1],"\n";
}
for ( $i = 0; $i <= ($pad*$sd); $i++ )
{
push @a1,$ax[0];
push @a2,$ax[@a-1];
}
push @a,@a1,@ax,@a2;
######## Gaussian Kernel###################################
my @x = map { 1 * $_ } (-1*($gw/2)) .. ($gw/2);
#### control the smoothness#######################
for ( $i = 0; $i < @x; $i++ )
{
$k[$i] = (1/sqrt(2*3.14159*$sd))* exp(-1*(($x[$i]*$x[$i])/(2*$sd*$sd)));
push @b,$k[$i];
}
######### Convolution noisy data with gaussian operator####
$la=@a; #length a
$lb=@b; #length b
$lab=$la+$lb-1; #length of result
for ( $i = 0+@x/2; $i < $lab-@x/2; $i++ )
{
$k=$i;
$y[$i]= 0;
for ( $j = 0; $j < $lb; $j++ ) #length b
{
if ($k>=0 && $k<$lab)
{
$y[$i] = $y[$i] + ($a[$k]*$b[$j]);
$k=$k-1;
}
}
push @yx, $y[$i];
push @yx_abs, abs($y[$i]);
}
$xxx2=@yx;
$scal=(average(@ax_abs))/(average(@yx_abs[$pad*$sd..$xxx2-($pad*$sd)-1]));
sub average { # mean of values in an array
my $sum = 0 ;
foreach $x (@_) {
$sum += $x ;
}
return $sum/@_ ;
}
open OUT2," > out.txt" or die "$!\n";
for ( $i = ($pad*$sd); $i < $xxx2-($pad*$sd); $i++ )
{
print OUT2 $scal*$yx[$i],"\n";
}
########## PLOT#############################################
open my $GP, '|-', 'gnuplot -persist';
print {$GP} <<'__GNUPLOT__';
plot 'in.txt' using 1 title 'BEFORE' with lines, 'out.txt' using 1 title 'AFTER' with lines
__GNUPLOT__
close $GP;
No comments:
Post a Comment