Binary Time constant 1.0 codec and decoder example

btc.py

#!/usr/bin/env python3

import array

BYTES = 2   # N bytes arthimetic
MAX = 2 ** (BYTES * 8 - 1) - 1
MIN = - (2 ** (BYTES * 8 - 1)) + 1

CHUNK= 1024

def predictive_btc(samples, soft = 21):
    """
    Encodes audio bytearray data with Binray Time constant codec to a BitStream
    in a list

    Keyword arguments:
        samples -- Signed 16 bit Audio data in a byte array
        soft - Softness constant of BTc. By default is 21

    """
    raw = array.array('h')
    raw.frombytes(samples)
    stream = []
    integrator = 0

    for sample in raw:
        # Reduces distance to max value in 1 / soft
        dist = (MAX - integrator) // soft
        highval = integrator + dist
        
        # Reduces distance to min value in 1 / soft
        dist = (integrator - MIN) // soft
        lowval = integrator - dist

        disthigh = abs(highval - sample)
        distlow = abs(lowval - sample)
        
        # Choose integrator with less diference to sample value
        if disthigh >= distlow:
            stream.append(0)
            integrator = lowval
        else:
            stream.append(1)
            integrator = highval
    
    return stream


def decode_btc(stream, soft = 21):
    """
    Decodes a Binary Time constant BitStream in Signed 16 bit Audio data in a
    ByteArray.

    Keywords arguments:
        stream -- List with the BitStream
        soft - Softness constant of BTc. By default is 21
    
    """

    audio = array.array('h')
    integrator = 0

    for bit in stream:
        if bit:     # Charge up
            integrator += (MAX - integrator) // soft
        else:       # Charge down
            integrator -= (integrator - MIN) // soft

        audio.append(int(integrator))

    return audio.tobytes()

# Main !
if __name__ == '__main__':
    try:
        import pyaudio
    except ImportError:
        print("Wops! We need PyAudio")
        sys.exit(0)

    import random
    import audioop
    import wave
    import sys
    import time
    from math import exp, log10

    random.seed()
    
    p = pyaudio.PyAudio() # Init PyAudio

    wf = wave.open(sys.argv[1], 'rb')
    print("Filename: %s" % sys.argv[1])
    Fm = wf.getframerate()
    print("Sample Rate: %d" % Fm)
    bits = wf.getsampwidth()
    channels = wf.getnchannels()

    samples = wf.readframes(wf.getnframes())    # Get all data from wave file
    wf.close()

    if bits != BYTES:   # Convert to Signed 16 bit data
        samples = audioop.lin2lin(samples, bits, BYTES)
        if bits == 1 and min(samples) >= 0:     # It was 8 bit unsigned !
            samples = audioop.bias(samples, BYTES, MIN)

    if channels > 1:    # Convert to Mono
        samples = audioop.tomono(samples, BYTES, 0.75, 0.25)
    
    # Normalize at 0.75
    maxsample = audioop.max(samples, BYTES)
    samples = audioop.mul(samples, BYTES, MAX * 0.75 / maxsample)

    soft = 21
    print("Softness value = %d" % soft)

    # Convert to Delta Modulation
    bitstream = predictive_btc(samples, soft)

    # Swap random bits (simulate bit errors)
    ERROR_RATE = 0 # 0.1
    BIT_PROB = 0.5
    tmp = max(BIT_PROB, 1- BIT_PROB)
    print("Error rate %f" % (ERROR_RATE * tmp))

    for i in range(len(bitstream)):
        if random.random() <= ERROR_RATE:
            if random.random() <= BIT_PROB:
                bitstream[i] = 0
            else:
                bitstream[i] = 1


    # Reconvert to PCM
    audio = decode_btc(bitstream, soft)

    # Play it!
    stream = p.open(format=p.get_format_from_width(BYTES), \
                        channels=1, rate=Fm, output=True)
    data = audio[:CHUNK]
    i = 0
    while i < len(audio):
        stream.write(data)
        i += CHUNK
        data = audio[i:min(i+CHUNK, len(audio))]

    time.sleep(0.5)
    
    stream.stop_stream()
    stream.close()

    p.terminate()

    s = 0.0
    n = 0.0
    for i in range(min(len(samples), len(audio))):
        s = s + samples[i]**2
        n = n +(samples[i] - audio[i])**2

    in_signal = float(s) / len(samples)
    ns_signal = float(n) / len(samples)

    snr = in_signal / ns_signal
    snr_db = 10 * log10(snr)
    print("SNR ratio %f (%f dB) " % (snr, snr_db))