a small example of the discrete cosine transform

dct.py

import math


def _dct(octave, t):
    return math.cos(math.pi * octave * t)


def _product(func, width, accum, octave):
    ''' compute the contribution of an octave '''
    p = 0.0
    for i in range(0, width):
        # find the time interpolant
        t = float(i) / float(width)
        # difference between real function and approximation
        d = func(t) - accum[i]
        # accumulate this correlation
        v = _dct(octave, t)
        p += v * d
    # correlation over sample count
    return p / float(width)


def deconstruct(func, width, num_coefs):
    accum = [0.0] * width
    coefs = [0.0] * num_coefs
    for i in range(0, num_coefs):
        coefs[i] = _product(func, width, accum, i)
        # add the discovered basis
        for j in range(0, width):
            accum[j] += _dct(i, float(j) / float(width)) * coefs[i]
    return coefs


def reconstruct(coefs, width):
    out = [0.0] * width
    for i, c in enumerate(coefs):
        for j in range(0, width):
            t = float(j) / float(width)
            out[j] += _dct(i, t) * c
    return out


def _square(t):
    ''' test is a square wave function '''
    return 1.0 if t >= 0.5 else -1.0


def _reduxsin(t):
    return float(int(math.sin(t * math.pi * 2.0) * 4.0)) / 4.0


def _write_file(samples):
    with file('out.csv', 'w') as fd:
        for i, v in enumerate(samples):
            fd.write('{}{:f}'.format(' ,' if i > 0 else '', v))


def _main():
    width = 128
    coefs = deconstruct(_reduxsin, width, 64)
    samples = reconstruct(coefs, width)
    _write_file(samples)


if __name__ == '__main__':
    _main()