perlでFFT

モジュールを使えば簡単



use Math::FFT;

my $PI = 3.1415926539;
my $N = 64;
my $in;

#set data
for (my $k=0; $k<$N; $k++) {
$in->[$k] = cos(2*$PI*$k/$N)+ sin(3*$PI*$k/$N);
}

#fft
my $fft = new Math::FFT($in);
my $coeff = $fft->rdft();
my $spectrum = $fft->spctrm;
my $original_data = $fft->invrdft($coeff);

#input data and inverse-fft
open FH, ">data.d";
for (my $k=0; $k<$N; $k++) {
$tmp1=0;
$tmp2=0;
for (my $j=1; $j<$N/2; $j++) {
$tmp1+=$coeff->[2*$j]*cos(2*$PI*$j*$k/$N);
$tmp2+=$coeff->[2*$j+1]*sin(2*$PI*$j*$k/$N);
}
$tmp=(($coeff->[0] +$coeff->[$N/2]*cos($PI*$k))/2 +$tmp1 +$tmp2)*2/$N;
print FH "$k\t $in->[$k] \t$tmp \t $original_data->[$k]\n";
}
close FH;

#coefficient
open FH, ">coeff.d";
for (my $k=1; $k<$N/2; $k++) {
print FH "$k\t $coeff->[2*$k] \t$coeff->[2*$k+1]\n";
}
close FH;

#power spectrum
open FH, ">spct.d";
for (my $k=1; $k<$N/4; $k++) {
my $tmp=($coeff->[2*$k]*$coeff->[2*$k]+ $coeff->[2*$k+1]*$coeff->[2*$k+1])/$N/$N*2;
print FH "$k\t $spectrum->[$k] \t$tmp\n";
}
close FH;


3/28追記
上記では，直流成分の係数は表示していないので注意