#!/usr/bin/env python
from __future__ import print_function, division

import argparse
import sys
try:
    from future_builtins import zip
except ImportError:     # either before Python2.6 or one of the Python 3.*
    try:
        from itertools import izip as zip
    except ImportError: # not an early Python, so must be Python 3.*
        pass

import numpy as np
import scipy
import scipy.signal

sampling_freq = 20      # if not read from input

times=[]
values=[]
for line in sys.stdin:
    if line.startswith('#'):
        print (line.rstrip())
        if line.startswith('# Recording every'):
           # extract from input: # Recording every 0.05 sec (20.0 Hz)
           fields = line.split('(')
           sampling_freq=float(fields[1].split()[0])
        continue
    fields = [x for x in map(float, line.split())]
    if len(fields)>=2:
        times.append(fields[0])
        values.append(fields[1])

# The band-pass filter will pass signals with frequencies between
# low_end_cutoff and high_end_cutoff
low_end_cutoff = 0.3    # Hz
lo_end_over_Nyquist = low_end_cutoff/(0.5*sampling_freq)
high_end_cutoff = 3.0 # Hz
hi_end_over_Nyquist = high_end_cutoff/(0.5*sampling_freq)

# If the bandpass filter gets ridiculously large output values (1E6 or more),
# the problem is numerical instability of the filter (probably from using a
# high sampling rate).  
# The problem can be addressed by reducing the order of the filter (first argument) from 5 to 2.
bess_b,bess_a = scipy.signal.iirfilter(5,
            Wn=[lo_end_over_Nyquist,hi_end_over_Nyquist],
            btype="bandpass", ftype='bessel')
bandpass = scipy.signal.filtfilt(bess_b,bess_a,values)

# The low-pass filter will pass signals with frequencies
# below low_end_cutoff
bess_b,bess_a = scipy.signal.iirfilter(5, Wn=[lo_end_over_Nyquist],
            btype="lowpass", ftype='bessel')
lowpass = scipy.signal.filtfilt(bess_b,bess_a,values)

print("# Bessel bandpass to {:.3g}Hz to {:.3g}Hz".format(low_end_cutoff, high_end_cutoff))
print("# Bessel lowpass {:.3g}Hz".format(low_end_cutoff))

for t,v,b,l in zip(times,values,bandpass,lowpass):
    print("{:.7f}\t{:.6f}\t{:.6f}\t{:.6f}".format(t,v,b,l))