from __future__ import division
from virgotools import *
from collections import defaultdict
from matplotlib.mlab import psd, window_hanning
import matplotlib.pyplot as plt


# get all calibration lines
vpm = VPM('DAQ')
freqs = defaultdict(list)

for proc in vpm.procs.itervalues():
    for fn, iline, line in parse_config_iter(proc.config):
        if line[0] == 'ACL_SINELINE_CH':
            freq = line[4]
        elif line[0] == 'ACL_SINEWAVE_CH':
            freq = line[6]
        else:
            continue
        chan = line[1]
        print proc.name, chan, freq
        
        # skip some uninteresting lines
        if proc.name == 'CALnoise' and 'perm' not in chan: continue
        if proc.name in ('ENVnoise', 'Acl_test'): continue
        if 'noise' in chan.lower(): continue
        if 'test' in chan.lower(): continue
        if 'MHz' in chan: continue  # derotation frequencies for dig demod 

        freqs[eval(freq)].append((proc.name, chan))  # use eval instead of float since there are fractions like 3/4
        
for k,v in sorted(freqs.iteritems()): 
    print k, v


# get data
chan = 'SDB2_B1s2_PD1_Blended'
gstart = 1177267897
dur = 600
data = getChannel(chan, gstart, dur)


# calculate ASD
twin = 40
nfft = int(data.fsample * twin)
PSD, freq = psd(data.data, NFFT=nfft, Fs=data.fsample, 
                window=window_hanning, noverlap=nfft/2)
ASD = PSD**0.5

def alias(fsample, f):
    return abs(f - fsample * round(f / fsample))


# make plot
label_y = 1e-2

plt.semilogy(freq, ASD)

for fline, sources in freqs.iteritems():
    if fline <= data.fsample / 2:
        color = 'r'
        label = 'normal line'
    else:
        color = 'green'
        fline = alias(data.fsample, fline)
        label = 'aliased line'
    plt.axvline(fline, color=color, linestyle='--', label=label)
        
    s = ', '.join('%s:%s' % source for source in sources)
    plt.text(fline, label_y, s,rotation=90, ha='left', va='bottom')

plt.xlabel('Frequency (Hz)')
plt.ylabel('%s (%s/rtHz)' % (chan, data.unit))
plt.title('Calibration lines in %s\nfrequency resolution %.2f Hz, using %i seconds of data at %s'
          % (chan, 1/twin, dur, gps2str(gstart)))
plt.show()