kiwifax.py

import array
import codecs
import logging
import os
import struct
import sys
import time
import math
import cmath
import traceback
from optparse import OptionParser

import kiwiclient
import png


# Known bugs and missing features:
# * No automatic LPM detection; useful when a station switches between 60 and 120
# * No automatic white level correction
# * When IQ input is used, the moving average would be beneficial to lower loise level


def dump_to_csv(filename, data, mode='a'):
    with open(filename, mode) as fp:
        for x in data:
            fp.write("%.6f," % x)
        fp.write("\n")


def clamp(x, xmin, xmax):
    if x < xmin:
        x = xmin
    if x > xmax:
        x = xmax
    return x

def norm_clamp(x, xmin, xmax):
    return (clamp(x, xmin, xmax) - xmin) / (xmax - xmin)

def fm_detect(X, prev, shift):
    vals = array.array('f')
    for x in X:
        y = shift + cmath.phase(x * prev.conjugate()) / math.pi
        vals.append(y)
        prev = x
    return vals


def dft_complex(input):
    width = len(input)
    output = []
    w1d = complex(0, -2 * math.pi / width)
    w1 = 0
    for k in xrange(width):
        X = 0
        w2d = cmath.exp(w1)
        w2 = complex(1, 0)
        for n in xrange(width):
            X += input[n] * w2
            w2 *= w2d
        output.append(X)
        w1 += w1d
    return output

def idft_complex(input):
    width = len(input)
    width_inv = 1.0 / width
    output = []
    w1d = complex(0, 2 * math.pi / width)
    w1 = 0
    for n in xrange(width):
        X = 0
        w2d = cmath.exp(w1)
        w2 = complex(1, 0)
        for k in xrange(width):
            X += input[k] * w2
            w2 *= w2d
        output.append(X * width_inv)
        w1 += w1d
    return output

def bitreverse_sort(input):
    output = list(input)

    half_length = len(input) // 2
    j = half_length
    for i in xrange(1, len(input) - 1):
        if i < j:
            t = output[j]
            output[j] = output[i]
            output[i] = t
        k = half_length
        while k <= j:
            j -= k
            k = k >> 1
        j += k

    return output

def log2(x):
    return math.frexp(x)[1] - 1

def fft_core(x):
    length = len(x)

    for l in xrange(1, log2(length) + 1):
        le = 1 << l
        le2 = le >> 1
        w = 2 * math.pi / le
        s = cmath.exp(complex(0, -w))
        u = complex(1, 0)
        for j in xrange(1, le2 + 1):
            for i in xrange(j - 1, length, le):
                o = i + le2
                t = x[o] * u
                x[o] = x[i] - t
                x[i] = x[i] + t
            u *= s

def fft_complex(input):
    x = bitreverse_sort(input)
    fft_core(x)
    return x

def ifft_complex(input):
    "Computes an inverse FFT transform for complex-valued input"
    x = bitreverse_sort(input)
    x = [ v.conjugate() for v in x ]
    fft_core(x)
    n_inv = 1.0 / len(x)
    x = [ v.conjugate() * n_inv for v in x ]
    return x

def power_db(input):
    nf = 1.0 / len(input)
    return [ 10 * math.log10(abs(x) * nf) for x in input ]

def peak_detect(data, thresh):
    data = array.array('f', data)
    peak_radius = 50
    peaks = []
    while True:
        peak_index = 0
        peak_value = data[peak_index]
        for i in xrange(1, len(data)):
            if peak_value < data[i]:
                peak_value = data[i]
                peak_index = i
        if peak_value < thresh:
            break
        peaks.append((peak_index, peak_value))
        for i in xrange(max(peak_index - peak_radius, 0), min(peak_index + peak_radius + 1, len(data))):
            data[i] = -999
    return peaks


class FMDetectorAtan2:
    def __init__(self):
        self._prev = complex(0)
    def process(self, samples):
        Y = array.array('f')
        prev = self._prev
        for x in samples:
            y = cmath.phase(x * prev.conjugate()) / math.pi
            Y.append(y)
            prev = x
        self._prev = prev
        return Y

class IQConverterDDC:
    """Convert audio samples to IQ: digital down-convert method"""
    def __init__(self, fc):
        "fc is the LO frequency divided by the sample rate"
        self._w = cmath.rect(1, -fc * 2 * math.pi)
        self._v = complex(1)
    def process(self, samples):
        Y = []
        for x in samples:
            Y.append(x * self._v)
            self._v *= self._w
        return Y

class IQConverterFFT:
    def __init__(self):
        pass
    def process(self, samples):
        X = fft_complex([ complex(x) for x in samples ])
        w = 1 + len(X) // 2
        Y = []
        for i in xrange(0, w):
            Y.append(X[i])
        for i in xrange(w, len(X)):
            Y.append(complex(1e-6))
        return ifft_complex(Y)


def interp_hermite(t, p0, p1, p2, p3):
    c0 = p1
    c1 = 0.5 * (p2 - p0)
    c2 = p0 - 2.5 * p1 + 2 * p2 - 0.5 * p3
    c3 = 0.5 * (p3 - p0) + 1.5 * (p1 - p2)
    return c0 + (t * (c1 + (t * (c2 + (t * c3)))))

class Interpolator:
    def __init__(self, factor):
        self._buffer = array.array('f')
        self._t = 0
        self.set_factor(factor)
    def set_factor(self, factor):
        self._dt = factor
    def extend(self, samples):
        for x in samples:
            self._buffer.append(float(x))
    def __iter__(self):
        return self
    def next(self):
        t_int = math.trunc(self._t)
        t_frac = self._t - t_int
        if t_int + 3 >= len(self._buffer):
            self._flush()
            raise StopIteration()
        self._t += self._dt
        return interp_hermite(t_frac, self._buffer[t_int], self._buffer[t_int + 1], self._buffer[t_int + 2], self._buffer[t_int + 3])
    def _flush(self):
        t_int = math.trunc(self._t)
        t_new = min(t_int, len(self._buffer))
        if t_new > 0:
            self._t -= t_new
            self._buffer = self._buffer[t_new:]

class FIRFilter:
    def __init__(self, kernel):
        self._kernel = kernel
        self._buffer = []
    def process(self, samples):
        self._buffer.extend(samples)
        Y = []
        i = 0
        while i + len(self._kernel) < len(self._buffer):
            y = 0
            for j in xrange(len(self._kernel)):
                y += self._buffer[i+j] * self._kernel[-j-1]
            Y.append(y)
            i += 1
        self._buffer = self._buffer[i:]
        return Y

def generate_sinc(fc, length):
    "Generates a sinc kernel"

    h = []
    w = 2 * math.pi * fc
    zf = (length - 1) / 2
    for i in xrange(0, length):
        x = i - zf
        if x == 0:
            h.append(w)
        else:
            h.append(math.sin(w * x) / x)
    return h

def generate_cosine_window_3(length, a, b, c, d):
    w = (2 * math.pi) / (length - 1)
    return [(a - b * math.cos(w * i) + c * math.cos(2 * w * i) -d * math.cos(3 * w * i)) for i in xrange(0, length)]

def generate_blackman_nuttall_window(length):
    return generate_cosine_window_3(length, 0.3635819, 0.4891775, 0.1365995, 0.0106411)

def apply_window(h, hw):
    if len(h) != len(hw):
        raise ValueError("vectors must have equal lengths")
    return [ h[i] * hw[i] for i in xrange(len(hw)) ]


def mapper_df_to_intensity(dfs, black_thresh, white_thresh):
    for x in dfs:
        yield norm_clamp(x, black_thresh, white_thresh)

class Histogram:
    def __init__(self, bins, xmin, xmax):
        self._min = xmin
        self._max = xmax
        self._bins = [ 0 for i in xrange(bins) ]
    def put(self, x):
        x = clamp(x, self._min, self._max)
        x = (x - self._min) / (self._max - self._min)
        i = int(x * (len(self._bins) - 1))
        self._bins[i] += 1
    def clear(self):
        for i in xrange(len(self._bins)):
            self._bins[i] = 0
    def get(self):
        s = 1.0 / sum(self._bins)
        return [ x * s for x in self._bins ]


# Let them have a name
RADIOFAX_WHITE_FREQ = 2300
RADIOFAX_BLACK_FREQ = 1500
RADIOFAX_STARTSTOP_FREQ = 1900
RADIOFAX_IOC576_START_TONE = 300
RADIOFAX_IOC288_START_TONE = 675
RADIOFAX_STOP_TONE = 450

class KiwiFax(kiwiclient.KiwiSDRSoundStream):
    def __init__(self, options):
        super(KiwiFax, self).__init__()
        self._options = options
        self._ioc = options.ioc
        self._lpm = options.lpm

        self._state = 'idle'

        self._use_iq = options.iq_stream

        self._iqconverter = None
        self._iqfir = None
        self._tuning_offset = options.force_offset
        self._ss_window_size = 4096
        self._startstop_buffer = []
        self._startstop_score = 0

        self._prevX = complex(0)
        self._phasing_count = 0
        self._resampler = None
        self._line_scale_factor = 1.0 - 1e-6 * options.sr_coeff
        self._rows = []
        self._pixel_buffer = array.array('f')
        # TODO: compute instead of hardcoding
        self._pixels_per_line = 1809
        # NOTE: Kyodo pages are ~8500px
        self._max_height = options.max_height

        self._histoa = Histogram(200, -0.1, +0.1)
        self._histob = Histogram(257, 0, 1)

        self._new_roll()
        if options.force:
            self._switch_state('printing')

    def _switch_state(self, new_state):
        logging.info("Switching to: %s", new_state)
        self._state = new_state
        if new_state == 'idle':
            self._startstop_score = 0
            self._noise_score = 0
            self._histoa.clear()
            self._histob.clear()
        elif new_state == 'starting':
            pass
        elif new_state == 'phasing':
            self._new_roll()
            self._phasing_count = 0
        elif new_state == 'printing':
            self._startstop_score = 0
        elif new_state == 'stopping':
            pass

    def _setup_rx_params(self):
        df = 1500
        if self._use_iq:
            # Tuned to the baseband
            bw = (RADIOFAX_WHITE_FREQ - RADIOFAX_BLACK_FREQ) / 2 + df
            self.set_mod('iq', -bw, +bw, self._options.frequency)
        else:
            # Tuned to USB (-1900 Hz)
            self.set_mod('usb', RADIOFAX_BLACK_FREQ - df, RADIOFAX_WHITE_FREQ + df, self._options.frequency - 1.9)
        # TODO: figure out proper AGC parameters
        self.set_agc(True)
        self.set_inactivity_timeout(0)
        self.set_name('kiwifax.py')
        # self.set_geo('Antarctica')

    def _on_sample_rate_change(self):
        sample_rate = float(self._sample_rate)
        # Precompute everything that depends on the SR
        self._iqconverter = IQConverterDDC(RADIOFAX_STARTSTOP_FREQ / sample_rate)
        filter_width = 450 # Hz
        filter_taps = 17
        self._iqfir = FIRFilter(apply_window(generate_sinc(filter_width / sample_rate, filter_taps), generate_blackman_nuttall_window(filter_taps)))
        # Start/stop detection params
        resolution = sample_rate / self._ss_window_size
        self._bin_size = resolution
        self._startstop_center_bin = self._ss_window_size // 2 + 0
        self._start576_delta = int(RADIOFAX_IOC576_START_TONE / resolution)
        self._start288_delta = int(RADIOFAX_IOC288_START_TONE / resolution)
        self._stop_delta = int(RADIOFAX_STOP_TONE / resolution)
        self._ss_width = int(0.5 * (RADIOFAX_WHITE_FREQ - RADIOFAX_STARTSTOP_FREQ) / resolution)
        logging.info("Start/stop center bin: %d; width: %d", self._startstop_center_bin, self._ss_width)
        logging.info("Start side bins: %d/%d; stop side bins: %d/%d",
            self._startstop_center_bin+self._start576_delta, self._startstop_center_bin-self._start576_delta,
            self._startstop_center_bin+self._stop_delta, self._startstop_center_bin-self._stop_delta)
        # NOTE: tone width is halved -- it should be precise anyway
        self._ss_tone_width = int(0.5 * 0.5 * (RADIOFAX_STOP_TONE - RADIOFAX_IOC576_START_TONE) / resolution)
        # Pixel output params
        samples_per_line = sample_rate * 60.0 / self._lpm
        resample_factor = (samples_per_line / self._pixels_per_line) * self._line_scale_factor
        self._resampler = Interpolator(resample_factor)
        logging.info("Resampling factor: %f", resample_factor)
        contrast = 0.01
        brightness = 0.02
        shift = 0.00
        self._white_level = (2 * (RADIOFAX_WHITE_FREQ - RADIOFAX_STARTSTOP_FREQ) / sample_rate) - contrast - brightness + shift
        self._black_level = (2 * (RADIOFAX_BLACK_FREQ - RADIOFAX_STARTSTOP_FREQ) / sample_rate) + contrast + shift
        self._fc_factor = 2 * self._bin_size / sample_rate

    def _process_audio_samples(self, seq, samples, rssi):
        k = 1 / 32768.0
        samples = [ x * k for x in samples ]
        samples = self._iqconverter.process(samples)
        self._process_samples(seq, samples, rssi)

    def _process_iq_samples(self, seq, samples, rssi, gps):
        k = 1 / 32768.0
        samples = [ x * k for x in samples ]
        self._process_samples(seq, samples, rssi)

    def _process_samples(self, seq, samples, rssi):
        logging.info('Block: %08x, RSSI: %04d %s', seq, rssi, self._state)
        if not self._use_iq:
            samples = self._iqfir.process(samples)
        self._process_startstop(samples)
        self._process_pixels(samples)

    def _startstop_score_update(self, updown):
        if updown:
            if self._startstop_score < 10:
                self._startstop_score += 1
        elif self._startstop_score > 0:
            self._startstop_score -= 2
            if self._startstop_score < 0:
                self._startstop_score = 0

    def _process_startstop(self, samples):
        self._startstop_buffer.extend(samples)
        # Snip out a window for start/stop processing
        # Window size defines the overall size of the window
        # Window shift defines how many samples are discarded after each iteration
        # This allows for overlapping FFTs thus increasing temporal resolution
        window_shift = self._ss_window_size / 2
        while len(self._startstop_buffer) >= self._ss_window_size:
            window = self._startstop_buffer[:self._ss_window_size]
            self._startstop_buffer = self._startstop_buffer[window_shift:]
            self._process_startstop_piece(window)

    def _process_startstop_piece(self, samples):
        # Compute the power spectrum
        samples = fft_complex(samples)
        P = power_db(samples)
        # DC "removal" for IQ
        if self._use_iq:
            P[0] = P[1]
        # Panoramize
        P1 = P[len(P)//2:]
        P1.extend(P[:len(P)//2])
        P = P1
        # DUMP POINT
        if self._options.dump_spectra and self._state != 'idle':
            dump_to_csv(self._output_name + '-ss.csv', P)
        # Assume noise floor is the median value + 5dB
        Px = P[2048-425:2048+425]
        Psorted = sorted(Px)
        nf_level = Psorted[len(Psorted) // 2] + 5.0
        pk_level = Psorted[-1]
        peaks = peak_detect(P, nf_level + 10)
        logging.info("Peaks: [%s]", ' '.join([ '%04d:%+05.1f' % (x[0], x[1]) for x in peaks ]))
        # For each peak, test if it's the one around the start/stop middle freq
        # For 4096-wide FFT: W=981 B=640 S=810 Start576=[682,939], Stop=[618,1002]
        detect_startstop = False
        detect_start576L = False
        detect_start576H = False
        detect_stopL = False
        detect_stopH = False
        # Classify the peaks
        for peak_bin, peak_power in peaks:
            # Try to classify the peak
            # Don't apply tuning correction for the start/stop center peak
            if math.fabs(peak_bin - self._startstop_center_bin) < self._ss_width:
                # NOTE: If force started, this doesn't get triggered properly
                if self._state in ('idle', 'starting'):
                    self._tuning_offset = self._startstop_center_bin - peak_bin
                detect_startstop = True
            else:
                peak_bin_relative = peak_bin + self._tuning_offset - self._startstop_center_bin
                if math.fabs(peak_bin_relative - self._stop_delta) < self._ss_tone_width:
                    detect_stopL = True
                if math.fabs(peak_bin_relative + self._stop_delta) < self._ss_tone_width:
                    detect_stopH = True
                if math.fabs(peak_bin_relative - self._start576_delta) < self._ss_tone_width:
                    detect_start576L = True
                if math.fabs(peak_bin_relative + self._start576_delta) < self._ss_tone_width:
                    detect_start576H = True
        detect_start576 = detect_startstop and detect_start576L and detect_start576H
        detect_stop = detect_startstop and detect_stopL and detect_stopH
        if self._state in ('idle', 'starting'):
            self._startstop_score_update(detect_start576)
        else:
            self._startstop_score_update(detect_stop)

        logging.info("NF=%05.1f PK=%05.1f  TO=%+04d/%+06.2fHz SS=%02d %s%s%s%s%s",
            nf_level, pk_level, self._tuning_offset, self._tuning_offset * self._bin_size,
            self._startstop_score,
            "sS"[detect_startstop], '-5'[detect_start576L], '-5'[detect_start576H], "xX"[detect_stopL], "xX"[detect_stopH])
        # Decide
        if self._state == 'idle':
            if self._startstop_score >= 10:
                logging.critical("START DETECTED")
                self._switch_state('starting')
        elif self._state == 'starting':
            if self._startstop_score < 3:
                self._switch_state('phasing')
        elif self._state == 'printing':
            if self._startstop_score >= 10:
                logging.critical("STOP DETECTED")
                self._switch_state('stopping')
        elif self._state == 'stopping':
            if self._startstop_score < 3:
                self._flush_rows()
                self._switch_state('idle')

    def _new_roll(self):
        self._rows = []
        ts = time.strftime('%Y%m%dT%H%MZ', time.gmtime())
        self._output_name = '%s_%d' % (ts, int(self._options.frequency * 1000))
        if self._options.station:
            self._output_name += '_' + self._options.station

    def _process_pixels(self, samples):
        if not self._state in ('phasing', 'printing', 'stopping'):
            return
        shift = self._tuning_offset * self._fc_factor
        pixels = fm_detect(samples, self._prevX, shift)
        self._prevX = samples[-1]
        # DUMP POINT
        if self._options.dump_pixels:
            dump_to_csv(self._output_name + '-px.csv', pixels)
        for x in pixels:
            self._histoa.put(x)
        # Remap the detected region into [0,1)
        pixels = array.array('f', mapper_df_to_intensity(pixels, self._black_level, self._white_level))
        for x in pixels:
            self._histob.put(x)
        # Scale and adjust pixel rate
        self._resampler.extend(pixels)
        self._pixel_buffer.extend(self._resampler)

        if self._state == 'phasing':
            self._process_phasing()
        else:
            # Cut into rows of pixels
            while len(self._pixel_buffer) >= self._pixels_per_line:
                row = self._pixel_buffer[:self._pixels_per_line]
                self._pixel_buffer = self._pixel_buffer[self._pixels_per_line:]
                self._process_row(row)

    def _process_phasing(self):
        # Count attempts at phasing to avoid getting stuck
        self._phasing_count += 1
        # Skip 3-4 lines; it seems phasing is not reliable when started right away
        if self._phasing_count <= 3:
            self._pixel_buffer = self._pixel_buffer[self._pixels_per_line:]
            return
        if self._phasing_count >= 100:
            logging.error("Phasing failed! Starting anyway")
            self._switch_state('printing')
            return
        # Do a moving average of the pixel intensity
        phasing_pulse_size = 90
        i = 0
        while i + phasing_pulse_size < len(self._pixel_buffer):
            s = 0
            for j in xrange(i, i + phasing_pulse_size):
                s += clamp(self._pixel_buffer[j], 0, 1)
            s /= phasing_pulse_size
            if s >= 0.85:
                self._pixel_buffer = self._pixel_buffer[i + phasing_pulse_size * 3 // 4:]
                logging.info("Phasing OK")
                self._switch_state('printing')
                break
            i += 1
        else:
            self._pixel_buffer = self._pixel_buffer[max(0, i - phasing_pulse_size):]

    def _process_row(self, row):
        pixels = array.array('B')
        for x in row:
            pixels.append(int(clamp(x, 0, 1) * 255))
        self._rows.append(pixels)
        if len(self._rows) % 16:
            return
        self._flush_rows()
        if len(self._rows) >= self._max_height:
            logging.info("Length exceeded; cutting the paper")
            self._switch_state('idle')

    def _flush_rows(self):
        if not self._rows:
            return
        while True:
            with open(self._output_name + '.png', 'wb') as fp:
                try:
                    png.Writer(len(self._rows[0]), len(self._rows), greyscale=True).write(fp, self._rows)
                    break
                except KeyboardInterrupt:
                    pass
        # DUMP POINT
        if self._options.dump_histo:
            dump_to_csv(self._output_name + '-hh.csv', self._histoa.get(), 'w')
            dump_to_csv(self._output_name + '-hh.csv', self._histob.get(), 'a')


KNOWN_CORRECTION_FACTORS = {
    'kiwisdr.northlandradio.nz:8073': {
        11030.00: -11.0,
    },
    'travelx.org:8073': { # +7.0
        7795.00: +3.0,
        9165.00: -5.0,
        13988.50: +3.0,
        16971.00: +4.0,
    },
    'travelx.org:8074': {
        7795.00: +0.0,
        9165.00: -14.0,
    },
    'reute.dyndns-remote.com:8073': {
        7880.00: -15.0,
        13882.50: -15.0,
    },
    'sarloutca.ddns.net:8073': {
        7880.00: -11.0,
        13882.50: -11.0,
    },
    'szsdr.ddns.net:8073': {
        9165.00: -11.0,
    },
    '72.130.191.200:8073': {
        9982.50: -13.0,
        11090.00: -13.0,
        16135.00: -13.0,
    }
}

def main():
    sys.stdout = codecs.getwriter('utf-8')(sys.stdout)

    parser = OptionParser()
    parser.add_option('-k', '--socket-timeout', '--socket_timeout',
                      dest='socket_timeout', type='int', default=10,
                      help='Timeout(sec) for sockets')
    parser.add_option('-s', '--server-host', '--server_host',
                      dest='server_host', type='string',
                      default='localhost', help='server host')
    parser.add_option('-p', '--server-port', '--server_port',
                      dest='server_port', type='int',
                      default=8073, help='server port (default 8073)')
    parser.add_option('-q', '--iq',
                      dest='iq_mode',
                      action='store_true', default=False,
                      help='IQ data mode')

    parser.add_option('-f', '--freq',
                      dest='frequency',
                      type='float', default=4610,
                      help='Frequency to tune to, in kHz (will be tuned down by 1.9kHz)')
    parser.add_option('--station', '--station',
                      dest='station',
                      type='string', default=None,
                      help='Station ID to be appended to file names')
    parser.add_option('-F', '--force-start',
                      dest='force',
                      action='store_true', default=False,
                      help='Force the decoding without waiting for start tone or phasing')
    parser.add_option('--force-offset', '--force_offset',
                      dest='force_offset',
                      type='int', default=0,
                      help='When force decoding, apply this tuning offset (bins).')
    parser.add_option('-i', '--ioc',
                      dest='ioc',
                      type='int', default=576,
                      help='Index of cooperation; default: 576.')
    parser.add_option('-l', '--lpm',
                      dest='lpm',
                      type='int', default=120,
                      help='Lines per minute; default: 120.')
    parser.add_option('--sr-coeff', '--sr_coeff',
                      dest='sr_coeff',
                      type='float', default=0,
                      help='Sample frequency correction, ppm; positive if the lines are too short; negative otherwise')
    parser.add_option('--max-height', '--max_height',
                      dest='max_height',
                      type='int', default=2300,
                      help='Maximum page height; default: 2300.')
    parser.add_option('--dump-spectra', '--dump-spectra',
                      dest='dump_spectra',
                      action='store_true', default=False,
                      help='Dump block spectra to a CSV file')
    parser.add_option('--dump-pixels', '--dump-pixels',
                      dest='dump_pixels',
                      action='store_true', default=False,
                      help='Dump row pixels to a CSV file')
    parser.add_option('--dump-histo', '--dump_histo',
                      dest='dump_histo',
                      action='store_true', default=False,
                      help='Dump pixel intensity histograms to a CSV file')
    parser.add_option('--iq-stream', '--iq_stream',
                      dest='iq_stream',
                      action='store_true', default=False,
                      help='EXPERIMENTAL: use IQ stream instead of audio')

    (options, unused_args) = parser.parse_args()

    # Setup logging
    fmtr = logging.Formatter('%(asctime)s %(levelname)s: %(message)s', '%Y%m%dT%H%MZ')
    fmtr.converter = time.gmtime
    fh = logging.FileHandler('log_%s_%d_%d.log' % (options.server_host, options.server_port, int(options.frequency * 1000)))
    fh.setLevel(logging.DEBUG)
    fh.setFormatter(fmtr)
    ch = logging.StreamHandler()
    ch.setLevel(logging.INFO)
    ch.setFormatter(fmtr)
    rootLogger = logging.getLogger()
    rootLogger.setLevel(logging.INFO)
    rootLogger.addHandler(fh)
    rootLogger.addHandler(ch)

    logging.critical('* * * * * * * *')
    logging.critical('Logging started')

    if options.sr_coeff == 0:
        server_identity = '%s:%d' % (options.server_host, options.server_port)
        try:
            coeffs = KNOWN_CORRECTION_FACTORS[server_identity]
            known_coeff = coeffs[options.frequency]
            options.sr_coeff = known_coeff
            logging.info('Applying known correction %f for host %s', known_coeff, server_identity)
        except KeyError:
            pass

    while True:
        recorder = KiwiFax(options)

        # Connect
        try:
            recorder.connect(options.server_host, options.server_port)
        except KeyboardInterrupt:
            break
        except Exception as e:
            traceback.print_exc()
            print "Failed to connect, sleeping and reconnecting"
            time.sleep(15)
            continue
        # Record
        try:
            recorder.run()
            break
        except (kiwiclient.KiwiTooBusyError, kiwiclient.KiwiBadPasswordError):
            print "Server too busy now, sleeping and reconnecting"
            time.sleep(15)
            continue
        except Exception as e:
            traceback.print_exc()
            break
    print "exiting"


if __name__ == '__main__':
    main()
# EOF