#!/usr/bin/env python

#      Filename: ftqviz.py
#        Author: Darren Hart <dvhltc@us.ibm.com>
#   Description: Plot the time and frequency domain plots of a times and
#                counts log file pair from the FTQ benchmark.
# Prerequisites: numpy, scipy, and pylab packages.  For debian/ubuntu:
#                o python-numeric
#                o python-scipy
#                o python-matplotlib
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# Copyright (C) IBM Corporation, 2007
#
# 2007-Aug-30:  Initial version by Darren Hart <dvhltc@us.ibm.com>


from numpy import *
from numpy.fft import *
from scipy import *
from pylab import *
from sys import *
from getopt import *

NS_PER_S  = 1000000000
NS_PER_MS = 1000000
NS_PER_US = 1000

def smooth(x, wlen):
    if x.size < wlen:
        raise ValueError, "Input vector needs to be bigger than window size."

    # reflect the signal to avoid transients... ?
    s = r_[2*x[0]-x[wlen:1:-1], x, 2*x[-1]-x[-1:-wlen:-1]]
    w = hamming(wlen)

    # generate the smoothed signal
    y = convolve(w/w.sum(), s, mode='same')
    # recenter the the smoothed signal over the originals (slide along x)
    y1 = y[wlen-1:-wlen+1]
    return y1


def my_fft(x, sample_hz):
    X = abs(fftshift(fft(x)))
    freq = fftshift(fftfreq(len(x), 1.0/sample_hz))
    return array([freq, abs(X)/len(x)])


def smooth_fft(timefile, countfile, sample_hz, wlen):
    # The higher the sample_hz, the larger the required wlen (used to generate
    # the hamming window).  It seems that each should be adjusted by roughly the
    # same factor
    ns_per_sample = NS_PER_S / sample_hz

    print "Interpolated Sample Rate: ", sample_hz, " HZ"
    print "Hamming Window Length: ", wlen

    t = fromfile(timefile, dtype=int64, sep='\n')
    x = fromfile(countfile, dtype=int64, sep='\n')

    # interpolate the data to achieve a uniform sample rate for use in the fft
    xi_len = (t[len(t)-1] - t[0])/ns_per_sample
    xi = zeros(xi_len)
    last_j = 0
    for i in range(0, len(t)-1):
        j = (t[i] - t[0])/ns_per_sample
        xi[j] = x[i]
        m = (xi[j]-xi[last_j])/(j-last_j)
        for k in range(last_j + 1, j):
            xi[k] = m * (k - last_j) + xi[last_j]
        last_j = j

    # smooth the signal (low pass filter)
    try:
        y = smooth(xi, wlen)
    except ValueError, e:
        exit(e)

    # generate the fft
    X = my_fft(xi, sample_hz)
    Y = my_fft(y, sample_hz)

    # plot the hamming window
    subplot(311)
    plot(hamming(wlen))
    axis([0,wlen-1,0,1.1])
    title(str(wlen)+" Point Hamming Window")
    
    # plot the signals
    subplot(312)
    ts = arange(0, len(xi), dtype=float)/sample_hz # time signal in units of seconds
    plot(ts, xi, alpha=0.2)
    plot(ts, y)
    legend(['interpolated', 'smoothed'])
    title("Counts (interpolated sample rate: "+str(sample_hz)+" HZ)")
    xlabel("Time (s)")
    ylabel("Units of Work")
    
    # plot the fft
    subplot(313)
    plot(X[0], X[1], ls='steps', alpha=0.2)
    plot(Y[0], Y[1], ls='steps')
    ylim(ymax=20)
    xlim(xmin=-3000, xmax=3000)
    legend(['interpolated', 'smoothed'])
    title("FFT")
    xlabel("Frequency")
    ylabel("Amplitude")

    show()


def usage():
        print "usage: "+argv[0]+" -t times-file -c counts-file [-s SAMPLING_HZ] [-w WINDOW_LEN] [-h]"


if __name__=='__main__':
    
    try:
        opts, args = getopt(argv[1:], "c:hs:t:w:")
    except GetoptError:
        usage()
        exit(2)
    
    sample_hz = 10000
    wlen = 25
    times_file = None
    counts_file = None
    for o, a in opts:
        if o == "-c":
            counts_file = a
        if o == "-h":
            usage()
            exit()
        if o == "-s":
            sample_hz = long(a)
        if o == "-t":
            times_file = a
        if o == "-w":
            wlen = int(a)

    if not times_file or not counts_file:
        usage()
        exit(1)

    smooth_fft(times_file, counts_file, sample_hz, wlen)