#======================================================================
#                    F F T . P L 
#                    doc: Fri Mar 12 09:20:33 1999
#                    dlm: Sun May 17 19:50:39 2015
#                    (c) 1999 A.M. Thurnherr
#                    uE-Info: 43 50 NIL 0 0 72 66 2 4 NIL ofnI
#======================================================================

# Notes:
#	It was found when rotary-analysing the FLAME current meters that
#	the sign of the frequencies returned was wrong. When investigating
#	the problem it was found that there are two conventions, one called
#	the engineering convention which is followed by Bendat & Piersol
#	[1971], Gonella [1972], and Mooers [1973] but not Numerical Recipes.
#	For real-valued functions it does not matter because the power
#	of the positive and negative frequencies are summed. The engineering
#	convention appears more sensible, however, because it actually leads
#	to the anticlockwise component to be reported as positive frequencies
#	which is consistent with exp(i phi) = cos phi + i sin phi and the
#	usual axes orientations.

# HISTORY:
#	Mar 12, 1999: - adapted from NR
#	Mar 13, 1999: - ``perlified'' (0-relative arrays; return value)
#	Mar 14, 1999: - cosmetic changes
#	Mar 15, 1999: - pad initial/final NaN values with 0es
#	Dec 08, 1999: - adapted for complex FFT
#	Dec 09, 1999: - continued
#	Dec 10, 1999: - BUG: < replaced by <= (no difference)
#	Dec 11, 1999: - investigated wrong frequency sign
#	Dec 14, 1999: - changed one-sided spectra to get half of the f=0 power
#	Mar 04, 2000: - require mean-removal when zero-padding is used
#				  - BUG (cosmetic)
#	Mar 05: 2000: - BUG: died even if no padding was done
#	Mar 08, 2000: - removed opt_r hard-coding
#	May 02, 2001: - BUG: RMS was really standard deviation estimator
#	Aug 29, 2003: - BUG: &icFFT did not calc sigma correctly whenever
#						 infield == outfield (e.g. [fftfilt] w/o -f)
#	Mar 31, 2004: - added cFFT_bufR()
#	Feb  9, 2005: - added &phase_pos(), &phase_neg()
#   Jul  1, 2006: - Version 3.3 [HISTORY]
#   Jul 24, 2006: - modified to use $PRACTICALLY_ZERO
#	May 15, 2015: - added $skip argument to cFFT()

#----------------------------------------------------------------------
# FOUR1 routine
#
# Notes:
#	- nan will abort FFT!
#	- added power of two assertions
#	- arrays are 0-relative
#	- isig should be set to -1 for FFT and to 1 for reverse FFT (see
#	  note above)
#
#----------------------------------------------------------------------

sub FOUR1($@) # ($isign, @data) => @fourier coefficients
{
	my($isign,@data) = @_;
	my($n,$mmax,$m,$j,$istep,$i);
	my($wtemp,$wr,$wpr,$wpi,$wi,$theta);
	my($tempr,$tempi);
	my($N) = @data / 2;

	$n = $N << 1;										# re-order input
	$j = 0;
	for ($i=0; $i<$n-1; $i+=2) {
		if ($j > $i) {
			$tempr = $data[$j];
			$tempi = $data[$j+1];
			$data[$j]   = $data[$i];
			$data[$j+1] = $data[$i+1];
			$data[$i]   = $tempr;
			$data[$i+1] = $tempi;
        }
        croak("$0 (fft.pl) $N is not a power of two\n")
        	if ($n % 2);
		$m = $n >> 1;
		while ($m >= 2 && $j >= $m) {
			$j -= $m;
			croak("$0 (fft.pl) $N is not a power of two\n")
	            if ($m % 2);
			$m >>= 1;
		}
		$j += $m;
	}

#	for ($i=0; $i<=$#data; $i+=2) {						# dump data
#		print(STDERR "$data[$i]$opt_O$data[$i+1]$opt_R")
#	}

	$mmax = 2;											# do FFT
	while ($n > $mmax) {
		$istep = $mmax << 1;
		$theta = $isign*(6.28318530717959/$mmax);
		$wtemp = sin(0.5*$theta);
		$wpr   = -2.0*$wtemp*$wtemp;
		$wpi   = sin($theta);
		$wr    = 1.0;
		$wi    = 0.0;
		for ($m=0; $m<$mmax-1; $m+=2) {
			for ($i=$m; $i<=$n-1; $i+=$istep) {
				$j     = $i + $mmax;
				$tempr = $wr*$data[$j]   - $wi*$data[$j+1];
				$tempi = $wr*$data[$j+1] + $wi*$data[$j];
				$data[$j]    = $data[$i]   - $tempr;
				$data[$j+1]  = $data[$i+1] - $tempi;
				$data[$i]   += $tempr;
				$data[$i+1] += $tempi;
			}
			$wtemp = $wr;
			$wr    = $wr*$wpr - $wi*$wpi    + $wr;
			$wi    = $wi*$wpr + $wtemp*$wpi + $wi;
		}
		$mmax = $istep;
	}

	return @data;
}

#----------------------------------------------------------------------
# TWOFFT routine
#
# Notes:
#	- arrays are 0-relative
#	- array references are used throughout; input arrays are 
#     in $ants_[][] style; ouput are simple arrays
#	- isign convention of NR is used here!!!
#----------------------------------------------------------------------

sub TWOFFT($$$$$$)
{
	my($data1R,$fnr1,$data2R,$fnr2,$fft1R,$fft2R) = @_;
	my($nn3,$nn2,$jj,$j);
	my($rep,$rem,$aip,$aim);
	my($n) = scalar(@{$data1R});

	$nn3 = 1 + ($nn2 = 2 + $n + $n);
	for ($j=1,$jj=2; $j<=$n; $j++,$jj+=2) {
		$fft1R->[$jj-2] = $data1R->[$j-1][$fnr1];
		$fft1R->[$jj-1] = $data2R->[$j-1][$fnr2];
	}
	@{$fft1R} = FOUR1(1,@{$fft1R});
	$fft2R->[0] = $fft1R->[1];
	$fft1R->[1] = $fft2R->[1] = 0;
	for ($j=3; $j<=$n+1; $j+=2) {
		$rep = 0.5 * ($fft1R->[$j-1] + $fft1R->[$nn2-$j-1]);
		$rem = 0.5 * ($fft1R->[$j-1] - $fft1R->[$nn2-$j-1]);
		$aip = 0.5 * ($fft1R->[$j] + $fft1R->[$nn3-$j-1]);
		$aim = 0.5 * ($fft1R->[$j] - $fft1R->[$nn3-$j-1]);
		$fft1R->[$j-1] = $rep;
		$fft1R->[$j]   = $aim;
		$fft1R->[$nn2-$j-1] =  $rep;
		$fft1R->[$nn3-$j-1] = -$aim;
		$fft2R->[$j-1] =  $aip;
		$fft2R->[$j]   = -$rem;
		$fft2R->[$nn2-$j-1] = $aip;
		$fft2R->[$nn3-$j-1] = $rem;
	}
}

#----------------------------------------------------------------------
# Interface to @ants_
#
# Notes:
#	- N (number of complex samples) calculated as next larger pwr-of-two
#	  if set to 0 on calling; otherwise set is used
#	- @ants_ padded with 0es to $N (deprecated in Hamming [1989])
#	- ditto initial and final missing values
#	- set ifnr to nan if input is purely real
#
#----------------------------------------------------------------------

sub cFFT($$$) { return cFFT_bufR(\@ants_,@_); }

sub cFFT_bufR($$$) # ($bufR, $rfnr, $ifnr, [$N[, $skip]]) => @coeff
{
	my($bufR,$fnr,$ifnr,$N,$skip) = @_;
	my(@data,$i,$r,$lastSet);

	unless ($N) {									# $N not set
		for ($N=1; $N <= $#ants_; $N <<= 1) {}		# next greater pwroftwo
		&antsInfo("(fft.pl) N set to $N")
			unless ($N == $#ants_+1);
	}
	for ($i=0,$r=$skip; $i<$N && $r<=$#ants_; $i++,$r++) {		# PAD
		last if (numberp($bufR->[$r][$fnr]) &&
				 (isnan($ifnr) || numberp($bufR->[$r][$ifnr])));
		$data[2*$i]   = 0;
		$data[2*$i+1] = 0;
	}
	$lastSet = $i - 1;
	&antsInfo("(fft.pl) WARNING: $i initial non-numbers padded with 0es!!!"),
		$padded=1 if ($i);
	while ($i<$N && $r<=$#ants_) {					# fill
		$i++,$r++,next
			unless (numberp($bufR->[$r][$fnr]) &&	# skip non-numbers 
				    (isnan($ifnr) || numberp($bufR->[$r][$ifnr])));
		croak("$0: (fft.pl) $lastSet, $i can't handle missing values ($bufR->[$lastSet+$skip+1][$fnr])!\n")
			if ($lastSet != $i-1);					# missing values
		$data[2*$i]   = $bufR->[$r][$fnr];			# real
		$data[2*$i+1] = isnan($ifnr) ? 0 : $bufR->[$r][$ifnr];	# imag
		$lastSet = $i;
		$i++; $r++;
	}
	&antsInfo("(fft.pl) WARNING: %d final non-numbers padded with 0es!!!",$i-$lastSet-1),
		$padded=1 if ($i > $lastSet+1);
	&antsInfo("(fft.pl) WARNING: padded with %d 0es to next pwr-of-two!!!",$N-$i),
		$padded=1 if ($i < $N);
	croak("$0: (fft.pl) refusing to pad with zeroes unless mean is removed (sorry)\n")
		if ($padded && !$FFT_ALLOW_ZERO_PADDING);
	$i = $lastSet + 1;
	while ($i < $N) {								# PAD
		$data[2*$i]   = 0;
		$data[2*$i+1] = 0;
		$i++;
	}
	return &FOUR1(-1,@data);
}

sub icFFT($$@) # ($ofnr, $tfnr, @coeff) => sigma
{
	my($ofnr,$tfnr,@coeff) = @_;
	my($N) = ($#coeff+1)/2;
	my(@val) = &FOUR1(1,@coeff);
	my($i);
	my($n) = 0;
	my($sumsq) = 0;
	my($mai) = 0;

	for ($i=0; $i<$N && $i<=$#ants_ ; $i++) {		# fill
		my($oldval) = $ants_[$i][$ofnr];
		push(@{$ants_[$i]},nan)						# pad empty fields
			while ($#{$ants_[$i]} < $tfnr);
		$ants_[$i][$tfnr] = $val[2*$i]/$N;			# real
		$mai = abs($val[2*$i+1])					# imag
			if (abs($val[2*$i+1]) > $mai);
		next unless (numberp($oldval));	# sigma
		$sumsq += ($oldval - $ants_[$i][$tfnr])**2;
		$n++;
	}
	&antsInfo("(fft.pl) WARNING: imaginary exponents (abs <= $mai) ignored")
		if ($mai > $PRACTICALLY_ZERO);
	&antsInfo("(fft.pl) WARNING: reducing data (%d -> %d)",$#ants_+1,$N)
		if ($i <= $#ants_);
	while ($i <= $#ants_) {
		$ants_[$i][$fnr] = nan;
		$i++;
    }
    return ($n > 1) ? sqrt($sumsq/($n-1)) : nan;
}

#----------------------------------------------------------------------
# Periodogram (p.421; (12.7.5) -- (12.7.6))
#----------------------------------------------------------------------

# Miscellaneous Notes:
#
# - there are N/2 + 1 values in the (unbinned) PSD (see NR)

# Notes regarding the effects of zero padding:
#
# - using zero-padding on a time series where the mean is not removed
#	can result in TOTALLY DIFFERENT RESULTS, try it e.g. with
#	temperatures!!! (A moment's thought will reveal why)
#
# - Because the total power is normalized to the mean squared amplitude
#	0-padded values depress the power; this was taken care of below by
#	normalizing the power by multiplying it with nrm=(nData+nZeroes)/nData;
#
# - this was checked using `avg -m' (in case of complex input the total
#	power is given by sqrt(Re**2+Im**2));
#
# - if zero-padding is used sqrt(nrm*P[0]) is mean value (done in [pgram])

# Notes on interpreting the power spectrum (Hamming [1989] & NR):
#
# - the frequency spectrum encompasses the frequencies between 0 (the
#	mean value) and the Nyquist frequency (1 / (2 x sampling interval))
#
# - higher frequencies are aliased into the power spectrum in a mirrored
#   way (e.g. noise tends to (linearly?) approach zero as f goes to inft;
#   the downsloping spectrum `hits' the Nyquist frequency, turns around
#   and continues falling towards the zero frequency, where it gets mirrored
# 	again => spectrum flattens towards Nyquist frequency
#
# - the sum over all P's (total power) is equal to the mean square value;
#	when one-sided spectra are used, P[0] and P[N/2] are counted doubly
#	and must be subtracted from the total; NB: the total power is reduced
#	if data are padded with 0es
#   
# - sqrt(P[0]) is an estimate for the mean value which is only accurate
#   if no zero-padding is perfomed; removing the mean will
#	strongly change the spectrum near the origin which might
#   or might not be a good thing, depending on the physics behind it (e.g.
#	it makes sense to remove the mean if a power spectrum from a temperature
#	record is calculated but not if flow velocity is used).
#   Removing higher order trends will also affect the spectrum but not
#   in such a simple fashion. Note that the problem is mainly restricted
#   to cases where the signal is in the low frequency and thus affected
#	by the strong changes in the spectrum there.

# Notes on the two-sided spectra:
# - the power of the mean flow (sqrt(P[0])) is non-rotary. To have the sum
#	of both one-sided spectra equal the two-sided one, each one-sided
#	spectrum gets half of the total value.
# - the same is true for the highest frequency; at the Nyquist frequency
#	every rotation is sampled exactly twice => polarization cannot be
#	determined (imagine a wheel with one spoke...)


sub pgram_onesided(@) # $nData,@C -> return @P
{
	my($nData,@C) = @_;
	my($N)    = ($#C+1) / 2;				# number of fourier comps
	my($Pfac) = $N**(-2) * $N/$nData;		# normalized to mean-sq amp
	my($k,@P);

	$P[0] = $Pfac * ($C[0]**2 + $C[1]**2);	# calc periodogram
	for ($k=1; $k<=$N/2-1; $k++) {
		$P[$k] = $Pfac * ($C[2*$k]**2 + $C[2*$k+1]**2 +
						  $C[2*($N-$k)]**2 + $C[2*($N-$k)+1]**2);
	}
	$P[$N/2] = $Pfac * ($C[2*($N/2)]**2 + $C[2*($N/2)+1]**2);
	return @P;
}

sub pgram_pos(@) # $nData,@C -> return @P
{
	my($nData,@C) = @_;
	my($N)    = ($#C+1) / 2;				# number of fourier comps
	my($Pfac) = $N**(-2) * $N/$nData;		# normalized to mean-sq amp
	my($k,@P);

	$P[0] = 0.5 * $Pfac * ($C[0]**2 + $C[1]**2);	# calc periodogram
	for ($k=1; $k<=$N/2-1; $k++) {
		$P[$k] = $Pfac * ($C[2*$k]**2 + $C[2*$k+1]**2);
	}
	$P[$N/2] = 0.5 * $Pfac * ($C[2*($N/2)]**2 + $C[2*($N/2)+1]**2);
	return @P;
}

sub pgram_neg(@) # $nData,@C -> return @P
{
	my($nData,@C) = @_;
	my($N)    = ($#C+1) / 2;				# number of fourier comps
	my($Pfac) = $N**(-2) * $N/$nData;		# normalized to mean-sq amp
	my($k,@P);

	$P[0] = 0.5 * $Pfac * ($C[0]**2 + $C[1]**2);	# calc periodogram
	for ($k=1; $k<=$N/2-1; $k++) {
		$P[$k] = $Pfac * ($C[2*($N-$k)]**2 + $C[2*($N-$k)+1]**2);
	}
	$P[$N/2] = 0.5 * $Pfac * ($C[2*($N/2)]**2 + $C[2*($N/2)+1]**2);
	return @P;
}

#------------------------------------------------------------------------
# Current Ellipses (Emery and Thomson, 5.6.4.2 (Rotary Component Spectra)
#------------------------------------------------------------------------

# Notes on phase (ellipsis inclination) calculations:
#	- comparing equations 5.6.42 and 5.6.43 in E+T shows a sign
#	  sign change of one of the terms (U_2k - V_1k in eqn 5.6.42a)
#	- for reasons of symmetry, this is unlikely to be correct
#	- in order to test this, I compared the M2 ellipses inclinations of
#	  a tidal analysis of BBTRE instruments #1 & #3 with the output
#	  and found that changing the sign in 5.6.43a brought the values
#	  into agreement (#1: tidal analysis 155+-3, pgram -25; #3: 136+-5,
#	  pgram -44)
#	- the changes are marked in the code by a (superfluous) + sign
#	- the sign change was later found to be consistent with the
#	  tidal analysis package manual by MGG Foreman (p.11)

sub phase_pos(@) # @C -> return @eps
{
	my(@C) = @_;
	my($N) = ($#C+1) / 2;				# number of fourier comps
	my($k,@eps);

	$eps[0] = 180 / $PI * atan2(+$C[1],$C[0]);
	for ($k=1; $k<=$N/2-1; $k++) {
		$eps[$k] = 180 / $PI * atan2(+$C[2*$k+1],$C[2*$k]);
	}
	$eps[$N/2] = 180 / $PI * atan2(+$C[2*($N/2)+1],$C[2*($N/2)]);
	return @eps;
}

sub phase_neg(@) # @C -> return @P
{
	my(@C) = @_;
	my($N) = ($#C+1) / 2;				# number of fourier comps
	my($k,@P);

	$eps[0] = 180 / $PI * atan2(+$C[1],$C[0]);
	for ($k=1; $k<=$N/2-1; $k++) {
		$eps[$k] = 180 / $PI * atan2(+$C[2*($N-$k)+1],$C[2*($N-$k)]);
	}
	$eps[$N/2] = 180 / $PI * atan2(+$C[2*($N/2)+1],$C[2*($N/2)]);
	return @eps;
}

#----------------------------------------------------------------------

1;