Descrete Fourier Transform in Perl

Sebuah sinyal dalam kawasan waktu  x[n]  di mana  n = 0,1,2,….,N-1 dan frekuensi X[k], dapat digambarkan dengan persamaan:

 

Teorema De Moivre’s menyatakan: 

  

Jika menggunakan Teorema De Moivre’s, maka persamaan di atas dapat dituliskan:

Komponen  Real:

Komponen Imajiner:

 

Gambar di bawah ini adalah sinyal dan spectrum frekuensi  yang diperoleh dengan menggunakan kode Perl terlampir:

#!/usr/bin/perl

#command: perl fft.pl input.txt sr
#sr in seconds

$file=$ARGV[0];   #get input file
$SR=$ARGV[1];     #sampling rate
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @x1, @row;
}
close( FILE );

open OUT1," > inputfft.txt" or die "$!\n";
for ($i=0; $i< @x1; $i++)
{
print OUT1 $i*$SR, "\t";
print OUT1 $x1[$i], "\n";
}

$n=@x1; #n input
$m=log($n)/log(2);
$k= int($m) + ($m != int($m)); #round up for power
$np=2**$k;   #n after padding

for ($i=0; $i< $np; $i++)
{
    push @x,1*$x1[$i];
}

$nx=@x;
$pi=3.14159;

for ( $k = 0; $k < $nx; $k++ )
{
#real part
$xre[$k]=0;
for ( $n = 0; $n < $nx; $n++ )
{
$xre[$k]+=$x[$n]*cos(((2*$pi*$k*$n)/$nx));
}

#imaginary part
$xim[$k]=0;
for ( $n = 0; $n < $nx; $n++ )
{
$xim[$k]-=$x[$n]*sin(((2*$pi*$k*$n)/$nx));
}
#power spectra
$P[$k]=sqrt((($xre[$k]*$xre[$k])+($xim[$k]*$xim[$k])));
push @Px, $P[$k];
}

open OUT," > outputfft.txt" or die "$!\n";
$maxP=max(@Px);
for ( $k = 0; $k <= (@Px/2); $k++ )
{
$f=($k/(@Px/2))*(1/(2*$SR));
print OUT $f, "\t";
print OUT $Px[$k]/$maxP, "\n";
}

sub max {
    splice(@_, ($_[0] > $_[1]) ? 1 : 0, 1);
    return ($#_ == 0) ? $_[0] : max(@_);
}

########## PLOT#############################################
open my $GP, '|-', 'gnuplot -persist';
print {$GP}  <<'__GNUPLOT__';

          set multiplot;      # get into multiplot mode
          set size 1,0.5; 
          set origin 0.0,0.5; 
    set xlabel 'Time (Second)'
    set ylabel 'Amplitude'
    set title 'TIME DOMAIN'
    set xtics font "Verdana,7"
    set ytics font "Verdana,7"
          plot 'inputfft.txt' using 1:2 title '' with lines;


          set origin 0.0,0.0; 
    set xlabel 'Frequency (Hz)'
    set ylabel 'Amplitude'
    set title 'FREQUENCY DOMAIN'
    set xtics font "Verdana,7"
    set ytics font "Verdana,7"
    plot 'outputfft.txt' using 1:2 title '' with lines
          unset multiplot    # exit multiplot mode

__GNUPLOT__

close $GP;