16acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak#!/usr/bin/env python 26acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 36acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# Filename: ftqviz.py 46acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# Author: Darren Hart <dvhltc@us.ibm.com> 56acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# Description: Plot the time and frequency domain plots of a times and 66acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# counts log file pair from the FTQ benchmark. 76acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# Prerequisites: numpy, scipy, and pylab packages. For debian/ubuntu: 86acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# o python-numeric 96acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# o python-scipy 106acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# o python-matplotlib 116acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# 126acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# This program is free software; you can redistribute it and/or modify 136acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# it under the terms of the GNU General Public License as published by 146acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# the Free Software Foundation; either version 2 of the License, or 156acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# (at your option) any later version. 166acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# 176acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# This program is distributed in the hope that it will be useful, 186acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# but WITHOUT ANY WARRANTY; without even the implied warranty of 196acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 206acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# GNU General Public License for more details. 216acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# 226acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# Copyright (C) IBM Corporation, 2007 236acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# 246acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak# 2007-Aug-30: Initial version by Darren Hart <dvhltc@us.ibm.com> 256acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 266acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 276acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakfrom numpy import * 286acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakfrom numpy.fft import * 296acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakfrom scipy import * 306acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakfrom pylab import * 316acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakfrom sys import * 326acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakfrom getopt import * 336acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 346acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakNS_PER_S = 1000000000 356acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakNS_PER_MS = 1000000 366acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakNS_PER_US = 1000 376acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 386acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakdef smooth(x, wlen): 396acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak if x.size < wlen: 406acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak raise ValueError, "Input vector needs to be bigger than window size." 416acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 426acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # reflect the signal to avoid transients... ? 436acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak s = r_[2*x[0]-x[wlen:1:-1], x, 2*x[-1]-x[-1:-wlen:-1]] 446acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak w = hamming(wlen) 456acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 466acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # generate the smoothed signal 476acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak y = convolve(w/w.sum(), s, mode='same') 486acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # recenter the the smoothed signal over the originals (slide along x) 496acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak y1 = y[wlen-1:-wlen+1] 506acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak return y1 516acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 526acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 536acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakdef my_fft(x, sample_hz): 546acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak X = abs(fftshift(fft(x))) 556acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak freq = fftshift(fftfreq(len(x), 1.0/sample_hz)) 566acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak return array([freq, abs(X)/len(x)]) 576acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 586acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 596acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakdef smooth_fft(timefile, countfile, sample_hz, wlen): 606acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # The higher the sample_hz, the larger the required wlen (used to generate 616acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # the hamming window). It seems that each should be adjusted by roughly the 626acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # same factor 636acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak ns_per_sample = NS_PER_S / sample_hz 646acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 656acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak print "Interpolated Sample Rate: ", sample_hz, " HZ" 666acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak print "Hamming Window Length: ", wlen 676acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 686acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak t = fromfile(timefile, dtype=int64, sep='\n') 696acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak x = fromfile(countfile, dtype=int64, sep='\n') 706acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 716acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # interpolate the data to achieve a uniform sample rate for use in the fft 726acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak xi_len = (t[len(t)-1] - t[0])/ns_per_sample 736acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak xi = zeros(xi_len) 746acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak last_j = 0 756acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak for i in range(0, len(t)-1): 766acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak j = (t[i] - t[0])/ns_per_sample 776acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak xi[j] = x[i] 786acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak m = (xi[j]-xi[last_j])/(j-last_j) 796acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak for k in range(last_j + 1, j): 806acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak xi[k] = m * (k - last_j) + xi[last_j] 816acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak last_j = j 826acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 836acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # smooth the signal (low pass filter) 846acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak try: 856acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak y = smooth(xi, wlen) 866acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak except ValueError, e: 876acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak exit(e) 886acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 896acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # generate the fft 906acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak X = my_fft(xi, sample_hz) 916acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak Y = my_fft(y, sample_hz) 926acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 936acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # plot the hamming window 946acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak subplot(311) 956acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak plot(hamming(wlen)) 966acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak axis([0,wlen-1,0,1.1]) 976acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak title(str(wlen)+" Point Hamming Window") 986acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 996acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # plot the signals 1006acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak subplot(312) 1016acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak ts = arange(0, len(xi), dtype=float)/sample_hz # time signal in units of seconds 1026acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak plot(ts, xi, alpha=0.2) 1036acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak plot(ts, y) 1046acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak legend(['interpolated', 'smoothed']) 1056acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak title("Counts (interpolated sample rate: "+str(sample_hz)+" HZ)") 1066acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak xlabel("Time (s)") 1076acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak ylabel("Units of Work") 1086acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1096acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak # plot the fft 1106acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak subplot(313) 1116acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak plot(X[0], X[1], ls='steps', alpha=0.2) 1126acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak plot(Y[0], Y[1], ls='steps') 1136acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak ylim(ymax=20) 1146acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak xlim(xmin=-3000, xmax=3000) 1156acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak legend(['interpolated', 'smoothed']) 1166acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak title("FFT") 1176acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak xlabel("Frequency") 1186acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak ylabel("Amplitude") 1196acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1206acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak show() 1216acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1226acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1236acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakdef usage(): 1246acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak print "usage: "+argv[0]+" -t times-file -c counts-file [-s SAMPLING_HZ] [-w WINDOW_LEN] [-h]" 1256acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1266acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1276acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modakif __name__=='__main__': 1286acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1296acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak try: 1306acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak opts, args = getopt(argv[1:], "c:hs:t:w:") 1316acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak except GetoptError: 1326acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak usage() 1336acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak exit(2) 1346acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1356acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak sample_hz = 10000 1366acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak wlen = 25 1376acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak times_file = None 1386acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak counts_file = None 1396acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak for o, a in opts: 1406acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak if o == "-c": 1416acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak counts_file = a 1426acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak if o == "-h": 1436acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak usage() 1446acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak exit() 1456acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak if o == "-s": 1466acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak sample_hz = long(a) 1476acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak if o == "-t": 1486acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak times_file = a 1496acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak if o == "-w": 1506acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak wlen = int(a) 1516acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1526acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak if not times_file or not counts_file: 1536acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak usage() 1546acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak exit(1) 1556acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak 1566acdc8efa73ceb0c3b515cd34c333d929e8b4273subrata_modak smooth_fft(times_file, counts_file, sample_hz, wlen) 157