#!/usr/bin/ruby

=begin
周波数解析を行う
=end

$: << "~/src/eclipse/autopilot/common/param/swig/build_SWIG"
$: << "~/src/eclipse/autopilot/common/ruby"

require 'FenrirMath'
require 'gnuplot'
require 'gnuplot_support'

$stderr.puts <<-__STRING__
Usage: #{__FILE__} input.csv [window_length = 256, freq_index = 1.. 20, skip = 16]
__STRING__

if ARGV.size < 1 then
  raise Exception::new("No file specified!!")
end

data = []
input_fname = ARGV.shift

def read_data(io)
  res = []
  io.each{|line|
    # 0列: k行目, 1列: 時刻, あとは入出力データ
    res << line.chop.split(/\s*,\s*/)[1..-1].collect{|item| item.to_f}
  }
  return res
end

if input_fname == '-' then
  data = read_data($stdin)
else
  open(input_fname){|io|
    data = read_data(io)
  }
end

WINDOW_LENGTH = eval(ARGV.shift) rescue 256
TARGET_SKIP = eval(ARGV.shift) rescue 16
TARGET_FREQ_INDEXES = eval(ARGV.shift) rescue 1..20

if data.size < WINDOW_LENGTH then
  raise Exception::new("Data too short!!")
end

top_index = 0
last_index = WINDOW_LENGTH - 1

fft_targets = []
(data[0].size - 1).times{|i|
  vec = FenrirMath::DVector::new
  vec.push(0)
  fft_targets << vec
}
data[0..(last_index-1)].each{|item|
  item[1..-1].each_with_index{|v, i|
    fft_targets[i].push(v)
  }
}

# http://solarwww.mtk.nao.ac.jp/sakamoto/memo/fft.htm

data_time = []  # X
data_freq = []  # Y
freq_step = 1.0 / (data[last_index][0] - data[top_index][0])
WINDOW_LENGTH.times{|i| data_freq << (freq_step * i)}
data_freq = data_freq[TARGET_FREQ_INDEXES] # 低周波しか興味がない
data_powers = [[]] * (data[0].size - 1) # Z

while last_index < data.size
  end_t, *values = data[last_index]
  values.each_with_index{|v, i|
    fft_targets[i].shift
    fft_targets[i].push(v)
  }
  start_t = data[top_index][0]
  if top_index % TARGET_SKIP == 0 then
    $stderr.puts "#{top_index}..#{last_index}"
    resolution_f = 1.0 / ((end_t - start_t) / WINDOW_LENGTH)
    fft_targets.each_with_index{|item, i|
      item = FenrirMath::WindowFunction::D_hamming(item)
      data_powers[i] << FenrirMath::FFT::D_fft(item)[TARGET_FREQ_INDEXES].collect{|complex|
        Math::log10(complex.abs2 / resolution_f) rescue 0.0
      }
    }
    data_time << (start_t + end_t) / 2
  end
  top_index += 1
  last_index += 1
end

eps_fname_base = File.basename(input_fname, ".*")
data_powers.each_with_index{|power, i|
  eps_fname = eps_fname_base + "_col_#{(?C + i).chr}.eps"
  Plotter::plot_basic_3d(eps_fname){|plot|
    
    #plot.xrange '[1:36]'
    #plot.yrange '[1.395:1.415]'
    #plot.set('zrange', '[0:2e+10]')
    
    plot.xlabel "'Time'"
    plot.ylabel "'Freq [Hz]'"
    #plot.set('zlabel', "'Power'")
    
    plot.set('pm3d map')
    #plot.set('palette', 'defined (400 "blue", 525 "white", 625 "red", 1200 "red")')
    #plot.set('dgrid3d', '41,71')
    #plot.set('hidden3d')
    
    #plot.set('ticslevel', 0.3)
    #plot.set('view 65,340')
    
    #p [data_time, data_freq, power]
    
    items = [
        Gnuplot::DataSet::new(DataSet3D::new([data_time, data_freq, power])){|ds|
          #ds.with = "dots" #"lines lw 2"
          #ds.title = "Analysis Result"
        }
      ]
      
    plot.data = items
  }
}