#!/usr/bin/env python

# (c) 2010 oswald berthold
#     during topology workshop, tm10

# - read usrp_spectrum_navigation log data, essentially the
#   output of usrp_spectrum_sense.py over a band of 800M - 2.46G
# - gps-tagged
# - multiple scans each with 2214*256 entries have to be split into
#   separate files
# - make separate plots for 1024 subbands chunks of the complete spectrum
#   for pdf / print compatibility
# - we just read one file

import time, sys
import numpy as np
import csv
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import matplotlib.ticker as ticker
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def dBm2mW(x):
    return 10**(x/10)/1000

#fname = "data/usrp_spectrum_navigation-20100205-141133-14.0-100.dat"
#fname = "data/usrp_spectrum_navigation-20100205-160603-2214.0-100.dat"


#reader = csv.reader(open(fname, "r"), delimiter="\t")
#rec = np.recfromcsv(fname, )
numfreq = 2214
fftsize = 256

chunklen = numfreq*fftsize
numchunk = 1
i = 1

freqs = np.zeros(chunklen)
spec_avg = np.zeros(chunklen)

# first, average over data
for i in range(1,32):
    fname = "../data/usn-05-16-%d.dat" % i
    print "fname:", fname
    f = open(fname, "r")
    for j in range(chunklen):
        line = f.readline()
        g = line.strip().split("\t", 5)
        #sys.stdout.write()
        #print g
        if i == 1: # once is enough
            freqs[j] = float(g[0])
        spec_avg[j] += float(g[1])

spec_avg /= 31.0 # normalise
ymax = dBm2mW(spec_avg.max())
#print spec_avg[:100]
print ymax

# then, make plots
#for i in range(chunklen / 1024):
numvalpp = 5000 # number of values per plot

for i in range(chunklen/numvalpp):
    
    # r = mlab.csv2rec(fname, delimiter="\t",
    #                  names=["freq", "dbm", "N", "E", "time", "alt", "nix"])

    # assert chunklen == r.size

    # # # 2D
    # #width = (2**15-1)/72 # 2214*256/72 # stupid max img size
    if i%4 == 0:
        width = 10
        fig = plt.figure(figsize=(width, 8))

    ax = fig.add_subplot(4,1,(i%4)+1)
    # ticks= np.arange(8e8, 2.6e9, 5e7) # 1e7
    # ax.set_xticks(ticks)
    ax.grid(True)
    ax.set_autoscaley_on(False)
    ax.set_ylim(0, ymax)
    #ax.set_xlabel("Freq [Hz]")
    ax.set_ylabel("{Power [mW]")

    if i%4 == 0:
        ax.set_title("Spectrum average (30 bins), 800M-2.46G")

    if i%4 == 3:
        ax.set_xlabel("Freq [Hz]")

    # # plt.plot(r['freq'][:256], \
    # #          r['dbm'][:256], axes=ax, clip_on=False)
    
    # # plt.plot(r['freq'][257:512], \
    # #          r['dbm'][257:512], axes=ax, clip_on=False)

    idxs = i*numvalpp # index-start
    idxe = (i+1)*numvalpp # index-end
    plt.plot(freqs[idxs:idxe], \
             dBm2mW(spec_avg[idxs:idxe]), \
             axes=ax, clip_on=False)

    # # plt.show()

    if i%4 == 3:
        fname2 = "../data/usn-05-16-spec-avg-%d-%d.png" % (idxs-(3*numvalpp), idxe)

        print "fname2: ", fname2
        plt.savefig(fname2, orientation="landscape", dpi=72)