# Make a bandpass filter to remove uninteresting frequencies
filt = sig.firwin(4*srate, [20.,480.], pass_zero=False, window='hann', fs=srate)
filt.resize(16*srate)

# Make the weighting filter - bandpass, weight by f^-7/6, and whiten
# This will give a divide by zero because of the f=0 point, ignore it
fweight = psd.frequencies**(-7./6.)*abs(rfft(filt))/sqrt(psd)
fweight[0] = 0. # Remove anything at f=0
norm = 1. / sqrt(sum(abs(fweight)**2) / (len(fweight)-1.))
fweight = norm * fweight.value

fwhiten = sqrt(2./srate)/sqrt(psd.value)
fwhiten[0] = 0. # Remove anything at f=0
td_whitener = sig.hann(16*srate)*roll(irfft(fwhiten), 8*srate)

full_filt = sig.hann(16*srate)*roll(irfft(fwhiten*fweight), 8*srate)

wdata = TimeSeries(sig.fftconvolve(data, full_filt, mode='same'),
                       sample_rate=data.sample_rate, epoch=data.epoch)[8*2048:-8*2048]

sum(abs(fweight)**2) / (16*2048/2.)

p1 = wdata.rms(8).plot()
p1.gca().set_ylim(0.8,1.2)

"""
Explanation: To get the normalization correct, it helps to break the filter into two pieces. The filters are based on the estimated (long) power spectrum denoted $S_{E}(f)$.
The first is the whitener, which is the inverse of the ASD in the frequency domain. The normalization that we want is that the whitened noise has variance one in the time domain.
$$\widetilde{W}(f) = \frac{f_{s}}{2} ~ \frac{1}{\sqrt{S_{E}(f)}}$$
The second is a filter that weights by the power in the whitened template. It looks like
$$\widetilde{H}(f) = \mathcal{N} ~ \frac{f^{-7/6}}{\sqrt{S_{E}(f)}}$$
The normalization is chosen so that the variance will be one if this filter is white noise which already has a variance of one. In the frequency domain, the sum of $|\widetilde{H}(f)|^2$ must equal the number of frequency bins in the filter (not including the one at $f = 0$).
End of explanation
"""


# Alternately, apply directly by FFT
# 1. Pad filter with zeros up to length of data
# 2. Multiply RFFT of data by abs(RFFT(filter))
# (Taking abs is the easy way to make sure a zero-phase filter has zero delay, even after padding)
# 3. Now take the IRFFT 

if False:
    wfilt_long = wfilt.copy()
    wfilt_long.resize(len(data))
    wdata = TimeSeries(irfft(rfft(data)*abs(rfft(wfilt_long))),
                       sample_rate=data.sample_rate, epoch=data.epoch)[8*2048:-8*2048]

# We can use the RMS to calculate the variance
# This is the same as integrating the spectrum, by Parseval's theorem
# But in practice, the spectrum may be better because we can use median or median-mean
# Note that I'm being lazy and not calculating the norm - just dividing by the median
# I'll do the norm next

my_rms = wdata.rms(8)
my_rms /= median(my_rms.value)

wdata.sample_rate

result = []
for idx in range(0, len(wdata)-8*2048, 4*2048):
    result.append(sqrt(median(wdata[idx:idx+8*2048].value**2)))
rmd = array(result)
rmd /= median(rmd)
rmd = TimeSeries(rmd, sample_rate=0.25, epoch=wdata.epoch)

p1 = (my_rms).plot()
p1.gca().plot(rmd)
p1.gca().set_ylim(0.9,2)
p1.gca().axvline(1170582798+20)

tconvert('2017-02-08 9:53:20')

"""
Explanation: The things below are not really needed. At this point, we've got a time series whose variance is the predicted variance of the SNR. Taking an RMS gives the predicted standard deviation.
End of explanation
"""


asd = data.asd(8,4,method='median')
weight = asd.frequencies**(-7./6.)/asd
weight[:8*20]=0.
weight[8*480:]=0.

sum(weight**2)/(sum(weight)**2)

sqrt(2*0.00040755575/(2))

sqrt(0.0001018889375)

asd2 = data.asd(1,0.5,method='median')
weight2 = asd2.frequencies**(-7./6.)/asd2
weight2[:1*20]=0.
weight2[1*480:]=0.

sum(weight2**2)/(sum(weight2)**2)

min(data.time

p1 = qtransform.plot()
p1.gca().plot()

1./(sum(weight**2)/(sum(weight)**2))

1./(sum(weight2**2)/(sum(weight2)**2))

"""
Explanation: Calculating the errors
This next section calculates the error on the PSD variation estimate, by modeling as a sum of weighted chi-squared variables. This can safely be ignored.
End of explanation
"""


