#===============================================
# MyFFT
#===============================================
package MyFFT;
use Exporter;
use Math::FFT;
@ISA = qw(Exporter Math::FFT);

#公開したいサブルーチン
#@EXPORT = qw(aa );

use strict;

sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;

	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub FFT
{
	my ($this) = @_;
	my $Direction = $this->{Direction};
	if($Direction =~ /inverse/i) {
		return $this->FFTInverse();
	}
	return $this->FFTForward();
}

sub FFTForward
{
	my ($this) = @_;

	my $pX    = $this->{pX};
	my $pData = $this->{pData};
	my $fft = new Math::FFT($pData);
	my ($n, $pCoeff);
	if($this->{Direction} =~ /real/i) {
		$n = @$pData;
		$pCoeff = $fft->rdft();
	}
	else {
		$n = int(@$pData / 2);
		$pCoeff = $fft->cdft();
	}
	$this->{pFowardCoeff} = $pCoeff;
#print "n=$n\n";

#print "x: $pX->[0], $pX->[1]\n";
	my $dx = $pX->[1] - $pX->[0];
	my $kmin = 1.0 / $dx / $n;
	my $kmax = 1.0 / $dx;
	my $kStep = ($kmax - $kmin) / ($n - 1);
	my (@k, @Areal, @Bimag);
	for(my $i = 0 ; $i < $n ; $i++) {
		if($this->{Direction} =~ /real/i) {
			last if($i >= $n / 2);
		}
		$k[$i] = $kmin + $i*$kStep;
		$Areal[$i] = $pCoeff->[2*$i];
		$Bimag[$i] = $pCoeff->[2*$i+1];
	}
	if($this->{Direction} =~ /real/i) {
		$Areal[$n] = $pCoeff->[$n/2];
		$Bimag[0] = 0.0;
	}
	$this->{pFowardK}     = \@k;
	$this->{pFowardAreal} = \@Areal;
	$this->{pFowardBimag} = \@Bimag;

	return (\@k, \@Areal, \@Bimag);
}

sub FFTInverse
{
	my ($this) = @_;

	my $pX    = $this->{pX};
	my $pData = $this->{pData};
	my $n = @$pData;
#print "n(inv)=$n\n";
	my $fft = new Math::FFT($pData);
	my $pCoeff;
	if($this->{Direction} =~ /real/i) {
		$pCoeff = $fft->invrdft($pData);
	}
	else {
		$pCoeff = $fft->invcdft($pData);
	}
	$this->{pFowardCoeff} = $pCoeff;

	my $dx = $pX->[1] - $pX->[0];
	my $kmin = 1.0 / $dx / $n;
	my $kmax = 1.0 / $dx;
	my $kStep = ($kmax - $kmin) / ($n - 1);
	my (@k, @Areal, @Breal);
	for(my $i = 0 ; $i < $n ; $i++) {
		$k[$i] = $kmin + $i*$kStep;
		if($this->{Direction} =~ /real/i) {
			$Areal[$i] = $pCoeff->[$i];
		}
		else {
			$Areal[$i] = $pCoeff->[2*$i];
			$Breal[$i] = $pCoeff->[2*$i+1];
		}
	}
	$this->{pInverseK}     = \@k;
	$this->{pInverseAreal} = \@Areal;
	$this->{pInverseAreal} = \@Breal;

	return (\@k, \@Areal, \@Breal);
}

sub SetData
{
# $Direction:    'forward',  'inverse', combined with 'complex', 'real', 'cosine' or 'sine'
# $AdjustMethod: 'truncate', 'fillzero'
	my $this = shift;
	my ($pX, $pReal, $pImag, $Direction, $AdjustMethod) = @_;
	if($Direction =~ /complex/i) {
		return $this->SetComplexData(@_);
	}
	return $this->SetRealData(@_);
}


sub SetComplexData
{
# $Direction:    'forward',  'inverse', combined with 'complex', 'real', 'cosine' or 'sine'
# $AdjustMethod: 'truncate', 'fillzero'
	my ($this, $pX, $pReal, $pImag, $Direction, $AdjustMethod) = @_;
	$this->{Direction}    = $Direction;
	$this->{AdjustMethod} = $AdjustMethod;

	my @x;
	my @data;
	my $nData = @$pReal;
#print "nData=$nData\n";
	for(my $i = 0 ; $i < $nData ; $i++) {
		last if($pX and (!defined $pX->[$i] or $pX->[$i] eq ''));
		if($pX) {
			$x[$i] = $pX->[$i];
		}
		else {
			$x[$i] = $i+1;
		}
		if($pReal) {
			$data[2*$i]   = $pReal->[$i];
		}
		else {
			$data[2*$i]   = 0.0;
		}
		if($pImag) {
			$data[2*$i+1] = $pImag->[$i];
		}
		else {
			$data[2*$i+1] = 0.0;
		}
	}
	$nData = int(@data / 2);
	my $lnn = int(log($nData)/log(2));
	my $n;
	if($AdjustMethod =~ /zero/i) {
		$n = 2**($lnn+1);
	}
	else {
		$n = 2**$lnn;
		my $n2 = 2*$n;
		@x    = @x[0..$n-1];
		@data = @data[0..$n2-1];
	}
#print "n=$n\n";

	$this->{pOriginalX}    = $pX;
	$this->{pOriginalReal} = $pReal;
	$this->{pOriginalImag} = $pImag;
	$this->{pX}            = \@x;
	$this->{pData}         = \@data;
	return ($n, \@x, \@data);
}

sub SetRealData
{
# $Direction:    'forward',  'inverse', combined with 'complex', 'real', 'cosine' or 'sine'
# $AdjustMethod: 'truncate', 'fillzero'
	my ($this, $pX, $pReal, $pImag, $Direction, $AdjustMethod) = @_;
	$this->{Direction}    = $Direction;
	$this->{AdjustMethod} = $AdjustMethod;

	my @x;
	my @data;
	my $nData = @$pReal;
#print "nData=$nData\n";
	for(my $i = 0 ; $i < $nData ; $i++) {
		last if($pX and (!defined $pX->[$i] or $pX->[$i] eq ''));

		if($pX) {
			$x[$i] = $pX->[$i];
		}
		else {
			$x[$i] = $i+1;
		}
		if($Direction =~ /forward/i) {
			if($pReal) {
				$data[$i]   = $pReal->[$i];
			}
			else {
				$data[$i]   = 0.0;
			}
		}
		else {
			if($pReal) {
				$data[2*$i]   = $pReal->[$i];
			}
			else {
				$data[2*$i]   = 0.0;
			}
			if($pImag) {
				$data[2*$i+1] = $pImag->[$i];
			}
			else {
				$data[2*$i+1] = 0.0;
			}
		}
	}
	if($Direction =~ /inverse/i) {
		$data[1] = $pReal->[$nData];
	}

	$nData = $nData;
	my $lnn = int(log($nData)/log(2));
	my $n;
	if($AdjustMethod =~ /zero/i) {
		$n = 2**($lnn+1);
	}
	else {
		$n = 2**$lnn;
		@x    = @x[0..$n-1];
		@data = @data[0..$n-1];
	}
#print "n=$n\n";

#	if($Direction =~ /inverse/i) {
#		$data[1] = 0.0;
#	}

	$this->{pOriginalX}    = $pX;
	$this->{pOriginalReal} = $pReal;
	$this->{pX}            = \@x;
	$this->{pData}         = \@data;
	return ($n, \@x, \@data);
}

1;