"""
apnea_sgram.py a11_Rtimes bphr expert cinc_path sgram

A hacked version of old record_view.py.  I created this
script to make a spectrogram figure for the book.
"""

import matplotlib, os, sys, re, scipy
a11_Rtimes, bphr, expert, cinc_path, sgram = sys.argv[1:6]
sys.path.append(cinc_path)
matplotlib.use('PDF')
import pylab, cinc2000
recordname = 'a11'

LOW    = 0.06
HIGH   = 0.14
W      = 7
expert = cinc2000.ReadIn(expert)

G_CM  = pylab.cm.hsv
#G_CM   = pylab.cm.gray

pylab.rc('figure', figsize=(5,3)) # same size as Gnuplot plots.
pylab.rc('lines', color='black')

# Read the data and do the calculations

bphr = file(bphr,'r')
HR_t = []
HR = []  # The heart rate
for line in bphr.readlines():
    part = line.split()
    HR_t.append(float(part[0]))
    HR.append(float(part[2]))
    # time from first column, low pass heart rate from third
Class_e = cinc2000.As_Ns_2_Nums(expert['a11'],val_N=-1,val_A=1)
Class_t = scipy.arange(0,len(Class_e))
#Fs = float(1)/6.0# samples per second
# Next get data for respiration spectrogram
data = []
for line in open(a11_Rtimes,'r').xreadlines():
    data.append(float(line)/100)
x = cinc2000.R_times2Dev(data,w=1) # Heart Rate Deviations
W = 6
NFFT=2**W
Fs = float(120)      # samples per minute
First = 4800-32      # 120*40
Last =  27000+40     # 120*225
Rxx, Rfreqs, Rbins, im = pylab.specgram(x[First:Last], NFFT=NFFT, Fs=Fs,
              window=pylab.window_hanning, noverlap=NFFT/2)
Rbins = Rbins + First/Fs
p = Rxx*Rxx
p = p.sum(axis=0)
p = 1/scipy.sqrt(p)
Rxx = Rxx*p

# Do the plotting
pylab.clf() # Clear figure window
ax1 = pylab.subplot(311)
pylab.plot(HR_t,HR)                   # Heart rate
pylab.ylim(-10,20)
yrng = scipy.arange(-5,20,5)
pylab.yticks(yrng, [ '$% 3.1f$' % l for l in yrng ])
pylab.xlim(40,225)
pylab.xticks(scipy.arange(50,250,50), visible=False)

# Plot respiration spectrogram
pylab.subplot(312,  sharex=ax1)
Z = -10*scipy.log10(Rxx[:-10,:]) #FixMe: The yaxis is wrong
Z =  scipy.flipud(Z)
extent = scipy.amin(Rbins), scipy.amax(Rbins), scipy.amin(Rfreqs),\
    scipy.amax(Rfreqs)/1.8
cmap = G_CM
Z = Z - Z.min()
vmin = LOW*Z.max()
vmax = HIGH*Z.max()
pylab.imshow(Z, cmap, extent=extent)
yrng = scipy.arange(0,35,10)
pylab.yticks(yrng, [ '$% 3.1f$' % l for l in yrng ])
pylab.xlim(40,225)
pylab.xticks(scipy.arange(50,250,50), visible=False)

pylab.subplot(313,  sharex=ax1)
pylab.plot(Class_t,Class_e)
pylab.ylim(-1.2,1.2)
pylab.xlim(40,225)
yrng = scipy.arange(-1,2,2)
pylab.yticks(yrng, ['$N$','$A$' ])
xrng = scipy.arange(50,250,50)
pylab.xticks(xrng,  [ '$%# 3i$' % l for l in xrng ])

pylab.savefig(sgram)

# Local Variables:
# mode: python
# End: