#!/usr/bin/perl -w  

#####################################################
# Program to convert astronomical spectra into waveforms, ready for conversion into
# sound. Uses the Fast Fourier Transform.
#
# Reads in an ASCII spectrum file, consisting of wavelengths and 
# fluxes (per unit frequency).
#
# Outputs a waveform file suitable for use with SoX.
#
# run as follows:
# ./specsound.pl <myspec.d > spec.dat
# ./sox spec.dat spec.wav
#
# Paul Francis, ANU, 7th Sept 2005.
#####################################################

use Math::FFT;
my $PI = 3.1415926539;

# Length of time (in 5.96 sec blocks) of sample clip. Integer
my $lclip = 2;

# FFT size. Must be 2**n.
my $N = 262144;
  
# Frequency of H-alpha (Hz): 261.63 = middle C.
my $freq_alpha = 261.63;    
  
# Set up the power spectrum. 
  
my ($series);
for (my $k=0; $k<$N; $k++) {
    $spec->[$k] = 0.0;
    $series->[2*$k] = 0.0;
    $series->[2*$k+1] = 0.0;
}


# Read in wavelength spectrum.

$nspec = 0;
while(<STDIN>){
   my @numbers = split;
   $wave[$nspec] = $numbers[0];
   $flux[$nspec] = abs($numbers[1]);
   $nspec++;
}


#Convert to frequency

for (my $k=0; $k < $nspec; $k++){
   $fnu[$nspec - $k - 1] = $flux[$k]*$wave[$k]*$wave[$k]/1.0E6;
   $dfreq[$nspec - $k - 1] = 5.96*$freq_alpha*6563.0/$wave[$k];
}


#Resample to frequency bin.

for (my $freq=0; $freq<$N; $freq++) {
   if ($freq > $dfreq[0] && $freq < $dfreq[$nspec-1]){
   
      $nd = 0;
      while($freq > $dfreq[$nd]){
        $nd++;
      }
      $spec->[$freq] += $fnu[$nd];
   }
}

# Convert from power spectrum to amplitude spectrum.

for (my $k=0; $k<$N; $k++) {
#    printf "%g %g \n", $k, $spec->[$k];
    $spec->[$k] = sqrt($spec->[$k]);
}

# Convert to complex number, adding random phases.

for (my $k=0; $k<$N; $k++) {
    my $phase = rand(2*$PI);
#    my $phase = $PI;
    $series->[2*$k] = cos($phase)*$spec->[$k];
    $series->[2*$k+1] = sin($phase)*$spec->[$k];
}

# Do the fourier transform

my $fft = new Math::FFT($series);
my $coeff = $fft->cdft();

# Normalise

my $max = 0;
my $min = 0;
my $norm = 1;
for (my $k=0; $k<$N; $k++) {
   if ($coeff->[2*$k] > $max){
      $max = $coeff->[2*$k];
   }
   if ($coeff->[2*$k] < $min){
      $min = $coeff->[2*$k];
   }
}
if ($max + $min > 0.0){
   $norm = $max;
} else {
   $norm = -1.0*$min;
}


# Output waveform.

print "; Sample Rate 44100\n"; # Correct value for garageband import.
my $nstart=0;
my $count = 0;
for (my $j=0; $j<$lclip; $j++){
   for (my $k=0; $k<$N; $k++) {
       my $temp = $coeff->[2*$k]/$norm;
       printf "%g %g \n", $count, $temp;    
       $count++;
   }
}