Saturday, January 24, 2015

Smoothing with Convolution


#!/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: