• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1#!/usr/bin/env python3
2
3#      Filename: ftqviz.py
4#        Author: Darren Hart <dvhltc@us.ibm.com>
5#   Description: Plot the time and frequency domain plots of a times and
6#                counts log file pair from the FTQ benchmark.
7# Prerequisites: numpy, scipy, and pylab packages.  For debian/ubuntu:
8#                o python-numeric
9#                o python-scipy
10#                o python-matplotlib
11#
12# This program is free software; you can redistribute it and/or modify
13# it under the terms of the GNU General Public License as published by
14# the Free Software Foundation; either version 2 of the License, or
15# (at your option) any later version.
16#
17# This program is distributed in the hope that it will be useful,
18# but WITHOUT ANY WARRANTY; without even the implied warranty of
19# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20# GNU General Public License for more details.
21#
22# Copyright (C) IBM Corporation, 2007
23#
24# 2007-Aug-30:  Initial version by Darren Hart <dvhltc@us.ibm.com>
25
26
27from numpy import *
28from numpy.fft import *
29from scipy import *
30from pylab import *
31from sys import *
32from getopt import *
33
34NS_PER_S  = 1000000000
35NS_PER_MS = 1000000
36NS_PER_US = 1000
37
38def smooth(x, wlen):
39    if x.size < wlen:
40        raise ValueError("Input vector needs to be bigger than window size.")
41
42    # reflect the signal to avoid transients... ?
43    s = r_[2*x[0]-x[wlen:1:-1], x, 2*x[-1]-x[-1:-wlen:-1]]
44    w = hamming(wlen)
45
46    # generate the smoothed signal
47    y = convolve(w/w.sum(), s, mode='same')
48    # recenter the smoothed signal over the originals (slide along x)
49    y1 = y[wlen-1:-wlen+1]
50    return y1
51
52
53def my_fft(x, sample_hz):
54    X = abs(fftshift(fft(x)))
55    freq = fftshift(fftfreq(len(x), 1.0/sample_hz))
56    return array([freq, abs(X)/len(x)])
57
58
59def smooth_fft(timefile, countfile, sample_hz, wlen):
60    # The higher the sample_hz, the larger the required wlen (used to generate
61    # the hamming window).  It seems that each should be adjusted by roughly the
62    # same factor
63    ns_per_sample = NS_PER_S / sample_hz
64
65    print("Interpolated Sample Rate: ", sample_hz, " HZ")
66    print("Hamming Window Length: ", wlen)
67
68    t = fromfile(timefile, dtype=int64, sep='\n')
69    x = fromfile(countfile, dtype=int64, sep='\n')
70
71    # interpolate the data to achieve a uniform sample rate for use in the fft
72    xi_len = (t[len(t)-1] - t[0])/ns_per_sample
73    xi = zeros(xi_len)
74    last_j = 0
75    for i in range(0, len(t)-1):
76        j = (t[i] - t[0])/ns_per_sample
77        xi[j] = x[i]
78        m = (xi[j]-xi[last_j])/(j-last_j)
79        for k in range(last_j + 1, j):
80            xi[k] = m * (k - last_j) + xi[last_j]
81        last_j = j
82
83    # smooth the signal (low pass filter)
84    try:
85        y = smooth(xi, wlen)
86    except ValueError as e:
87        exit(e)
88
89    # generate the fft
90    X = my_fft(xi, sample_hz)
91    Y = my_fft(y, sample_hz)
92
93    # plot the hamming window
94    subplot(311)
95    plot(hamming(wlen))
96    axis([0,wlen-1,0,1.1])
97    title(str(wlen)+" Point Hamming Window")
98
99    # plot the signals
100    subplot(312)
101    ts = arange(0, len(xi), dtype=float)/sample_hz # time signal in units of seconds
102    plot(ts, xi, alpha=0.2)
103    plot(ts, y)
104    legend(['interpolated', 'smoothed'])
105    title("Counts (interpolated sample rate: "+str(sample_hz)+" HZ)")
106    xlabel("Time (s)")
107    ylabel("Units of Work")
108
109    # plot the fft
110    subplot(313)
111    plot(X[0], X[1], ls='steps', alpha=0.2)
112    plot(Y[0], Y[1], ls='steps')
113    ylim(ymax=20)
114    xlim(xmin=-3000, xmax=3000)
115    legend(['interpolated', 'smoothed'])
116    title("FFT")
117    xlabel("Frequency")
118    ylabel("Amplitude")
119
120    show()
121
122
123def usage():
124        print("usage: "+argv[0]+" -t times-file -c counts-file [-s SAMPLING_HZ] [-w WINDOW_LEN] [-h]")
125
126
127if __name__=='__main__':
128
129    try:
130        opts, args = getopt(argv[1:], "c:hs:t:w:")
131    except GetoptError:
132        usage()
133        exit(2)
134
135    sample_hz = 10000
136    wlen = 25
137    times_file = None
138    counts_file = None
139    for o, a in opts:
140        if o == "-c":
141            counts_file = a
142        if o == "-h":
143            usage()
144            exit()
145        if o == "-s":
146            sample_hz = int(a)
147        if o == "-t":
148            times_file = a
149        if o == "-w":
150            wlen = int(a)
151
152    if not times_file or not counts_file:
153        usage()
154        exit(1)
155
156    smooth_fft(times_file, counts_file, sample_hz, wlen)
157