#!/usr/bin/ruby -Ks

require 'complex'

# FFT
class FFT

	def FFT.fft(data)

		#まず何回2で割り切れるか数える
		power2 = 0
		size = data.size
		while (size & 0b1) == 0
			power2 += 1
			size >>= 1
		end
		
		# FFTをする前に並べ替え
		border = data.size
		sorted = Array::new([data])
		(0...power2).each{|i|
			border >>= 1
			temp_plus = []
			temp_minus = []
			sorted.each{|element|
				sorted_plus = []
				sorted_minus = []
				(0...border).each{|j|
					sorted_plus << (element[j] + element[j + border])
					theta = 2 * j * Math::PI / element.size
					sorted_minus << (element[j] - element[j + border]) * Complex(Math::cos(theta), -Math::sin(theta))
				}
				temp_plus << sorted_plus
				temp_minus << sorted_minus
			}
			sorted = temp_plus + temp_minus
			#p sorted
		}

		# FFTする
		result = Array::new()
		base = 0
		(0...(data.size / sorted.size)).each{|j|
			k = j * (1 << power2)
			sorted.each_with_index{|element, i|
				result[base + i] = ft(element, k)
				#puts "i: #{i}; j: #{j}; k: #{k}"
			}
			base += (1 << (power2 * 2));
	        base %= data.size();
		}

		return result
	end

	# X_{k} = \sum_{n=0}^{N-1} x_{n} \exp^{\frac{2 \pi n k}{N}}
	def FFT.ft(data, k)
		result = 0
		data.each_with_index{|v, n|
			theta = 2 * Math::PI * n * k / data.size
			result += v * Complex::new(Math::cos(theta), -Math::sin(theta))
		}
		return result
	end
end

if __FILE__ == $0 then

DATA_FILE = './data.csv'	# データファイル名
DATA_FREQ = 5120.0			# データ取り込み周波数

#ファイルを開いてデータを取り込む
data = []
open(DATA_FILE){|src|
	src.each{|line|
		data << line.chop.to_f
	}
}

=begin
data = [1,2,1,2,1,2,1,2]
(0...data.size).each{|j|
	hoge = 0
	(0...data.size).each{|i|
		hoge += data[i] * (Math::E ** Complex::new(0, - 2 * Math::PI * i * j / data.size ))
	}
	p hoge
}
=end

step_freq = DATA_FREQ / data.size
rrs2 = 0
puts "Freq(Hz),Real,Imaginary,PSD,RRS"
FFT::fft(data).each_with_index{|c, i|
	freq = step_freq * i
	# 大域を50Hzから2KHzに制限
	if (50..2000).include?(freq) then
		psd = c.abs2 * 2 / data.size * (1 / DATA_FREQ)
		puts "#{freq},#{c.real},#{c.image},#{psd},#{Math::sqrt(rrs2 += psd * step_freq)}"
	end
}

end