transforms.py

vlc stand 0.1.99 ist sehr einfach

# Lossy compression algorithms very often make use of DCT or DFT calculations,
# or variations of these calculations. This file is intended to be a short
# reference about classical DCT and DFT algorithms.


import math
import cmath

pi = math.pi
sin = math.sin
cos = math.cos
sqrt = math.sqrt

def exp_j (alpha):
    return cmath.exp (alpha * 1j)

def conjugate (c):
    c = c + 0j
    return c.real - 1j * c.imag

def vector (N):
    return [0j] * N


# Let us start withthe canonical definition of the unscaled DFT algorithm :
# (I can not draw sigmas in a text file so I'll use python code instead)  :)

def W (k, N):
    return exp_j ((-2*pi*k)/N)

def unscaled_DFT (N, input, output):
    for o in range(N):		# o is output index
	output[o] = 0
	for i in range(N):
	    output[o] = output[o] + input[i] * W (i*o, N)

# This algorithm takes complex input and output. There are N*N complex
# multiplications and N*(N-1) complex additions.


# Of course this algorithm is an extremely naive implementation and there are
# some ways to use the trigonometric properties of the coefficients to find
# some decompositions that can accelerate the calculation by several orders
# of magnitude... This is a well known and studied problem. One of the
# available explanations of this process is at this url :
# www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html


# Let's start with the radix-2 decimation-in-time algorithm :

def unscaled_DFT_radix2_time (N, input, output):
    even_input = vector(N/2)
    odd_input = vector(N/2)
    even_output = vector(N/2)
    odd_output = vector(N/2)

    for i in range(N/2):
	even_input[i] = input[2*i]
	odd_input[i] = input[2*i+1]

    unscaled_DFT (N/2, even_input, even_output)
    unscaled_DFT (N/2, odd_input, odd_output)

    for i in range(N/2):
	odd_output[i] = odd_output[i] * W (i, N)

    for i in range(N/2):
	output[i] = even_output[i] + odd_output[i]
	output[i+N/2] = even_output[i] - odd_output[i]

# This algorithm takes complex input and output.

# We divide the DFT calculation into 2 DFT calculations of size N/2
# We then do N/2 complex multiplies followed by N complex additions.
# (actually W(0,N) = 1 and W(N/4,N) = -j so we can skip a few of these complex
# multiplies... we will skip 1 for i=0 and 1 for i=N/4. Also for i=N/8 and for
# i=3*N/8 the W(i,N) values can be special-cased to implement the complex
# multiplication using only 2 real additions and 2 real multiplies)

# Also note that all the basic stages of this DFT algorithm are easily
# reversible, so we can calculate the IDFT with the same complexity.


# A varient of this is the radix-2 decimation-in-frequency algorithm :

def unscaled_DFT_radix2_freq (N, input, output):
    even_input = vector(N/2)
    odd_input = vector(N/2)
    even_output = vector(N/2)
    odd_output = vector(N/2)

    for i in range(N/2):
	even_input[i] = input[i] + input[i+N/2]
	odd_input[i] = input[i] - input[i+N/2]

    for i in range(N/2):
	odd_input[i] = odd_input[i] * W (i, N)

    unscaled_DFT (N/2, even_input, even_output)
    unscaled_DFT (N/2, odd_input, odd_output)

    for i in range(N/2):
	output[2*i] = even_output[i]
	output[2*i+1] = odd_output[i]

# Note that the decimation-in-time and the decimation-in-frequency varients
# have exactly the same complexity, they only do the operations in a different
# order.

# Actually, if you look at the decimation-in-time varient of the DFT, and
# reverse it to calculate an IDFT, you get something that is extremely close
# to the decimation-in-frequency DFT algorithm...


# The radix-4 algorithms are slightly more efficient : they take into account
# the fact that with complex numbers, multiplications by j and -j are also
# "free"... i.e. when you code them using real numbers, they translate into
# a few data moves but no real operation.

# Let's start with the radix-4 decimation-in-time algorithm :

def unscaled_DFT_radix4_time (N, input, output):
    input_0 = vector(N/4)
    input_1 = vector(N/4)
    input_2 = vector(N/4)
    input_3 = vector(N/4)
    output_0 = vector(N/4)
    output_1 = vector(N/4)
    output_2 = vector(N/4)
    output_3 = vector(N/4)
    tmp_0 = vector(N/4)
    tmp_1 = vector(N/4)
    tmp_2 = vector(N/4)
    tmp_3 = vector(N/4)

    for i in range(N/4):
	input_0[i] = input[4*i]
	input_1[i] = input[4*i+1]
	input_2[i] = input[4*i+2]
	input_3[i] = input[4*i+3]

    unscaled_DFT (N/4, input_0, output_0)
    unscaled_DFT (N/4, input_1, output_1)
    unscaled_DFT (N/4, input_2, output_2)
    unscaled_DFT (N/4, input_3, output_3)

    for i in range(N/4):
	output_1[i] = output_1[i] * W (i, N)
	output_2[i] = output_2[i] * W (2*i, N)
	output_3[i] = output_3[i] * W (3*i, N)

    for i in range(N/4):
	tmp_0[i] = output_0[i] + output_2[i]
	tmp_1[i] = output_0[i] - output_2[i]
	tmp_2[i] = output_1[i] + output_3[i]
	tmp_3[i] = output_1[i] - output_3[i]

    for i in range(N/4):
	output[i] = tmp_0[i] + tmp_2[i]
	output[i+N/4] = tmp_1[i] - 1j * tmp_3[i]
	output[i+N/2] = tmp_0[i] - tmp_2[i]
	output[i+3*N/4] = tmp_1[i] + 1j * tmp_3[i]

# This algorithm takes complex input and output.

# We divide the DFT calculation into 4 DFT calculations of size N/4
# We then do 3*N/4 complex multiplies followed by 2*N complex additions.
# (actually W(0,N) = 1 and W(N/4,N) = -j so we can skip a few of these complex
# multiplies... we will skip 3 for i=0 and 1 for i=N/8. Also for i=N/8
# the remaining W(i,N) and W(3*i,N) multiplies can be implemented using only
# 2 real additions and 2 real multiplies. For i=N/16 and i=3*N/16 we can also
# optimise the W(2*i/N) multiply this way.

# If we wanted to do the same decomposition with one radix-2 decomposition
# of size N and 2 radix-2 decompositions of size N/2, the total cost would be
# N complex multiplies and 2*N complex additions. Thus we see that the
# decomposition of one DFT calculation of size N into 4 calculations of size
# N/4 using the radix-4 algorithm instead of the radix-2 algorithm saved N/4
# complex multiplies...


# The radix-4 decimation-in-frequency algorithm is similar :

def unscaled_DFT_radix4_freq (N, input, output):
    input_0 = vector(N/4)
    input_1 = vector(N/4)
    input_2 = vector(N/4)
    input_3 = vector(N/4)
    output_0 = vector(N/4)
    output_1 = vector(N/4)
    output_2 = vector(N/4)
    output_3 = vector(N/4)
    tmp_0 = vector(N/4)
    tmp_1 = vector(N/4)
    tmp_2 = vector(N/4)
    tmp_3 = vector(N/4)

    for i in range(N/4):
	tmp_0[i] = input[i] + input[i+N/2]
	tmp_1[i] = input[i+N/4] + input[i+3*N/4]
	tmp_2[i] = input[i] - input[i+N/2]
	tmp_3[i] = input[i+N/4] - input[i+3*N/4]

    for i in range(N/4):
	input_0[i] = tmp_0[i] + tmp_1[i]
	input_1[i] = tmp_2[i] - 1j * tmp_3[i]
	input_2[i] = tmp_0[i] - tmp_1[i]
	input_3[i] = tmp_2[i] + 1j * tmp_3[i]

    for i in range(N/4):
	input_1[i] = input_1[i] * W (i, N)
	input_2[i] = input_2[i] * W (2*i, N)
	input_3[i] = input_3[i] * W (3*i, N)

    unscaled_DFT (N/4, input_0, output_0)
    unscaled_DFT (N/4, input_1, output_1)
    unscaled_DFT (N/4, input_2, output_2)
    unscaled_DFT (N/4, input_3, output_3)

    for i in range(N/4):
	output[4*i] = output_0[i]
	output[4*i+1] = output_1[i]
	output[4*i+2] = output_2[i]
	output[4*i+3] = output_3[i]

# Once again the complexity is exactly the same as for the radix-4
# decimation-in-time DFT algorithm, only the order of the operations is
# different.


# Now let us reorder the radix-4 algorithms in a different way :

#def unscaled_DFT_radix4_time (N, input, output):
#   input_0 = vector(N/4)
#   input_1 = vector(N/4)
#   input_2 = vector(N/4)
#   input_3 = vector(N/4)
#   output_0 = vector(N/4)
#   output_1 = vector(N/4)
#   output_2 = vector(N/4)
#   output_3 = vector(N/4)
#   tmp_0 = vector(N/4)
#   tmp_1 = vector(N/4)
#   tmp_2 = vector(N/4)
#   tmp_3 = vector(N/4)
#
#   for i in range(N/4):
#	input_0[i] = input[4*i]
#	input_2[i] = input[4*i+2]
#
#   unscaled_DFT (N/4, input_0, output_0)
#   unscaled_DFT (N/4, input_2, output_2)
#
#   for i in range(N/4):
#	output_2[i] = output_2[i] * W (2*i, N)
#
#   for i in range(N/4):
#	tmp_0[i] = output_0[i] + output_2[i]
#	tmp_1[i] = output_0[i] - output_2[i]
#
#   for i in range(N/4):
#	input_1[i] = input[4*i+1]
#	input_3[i] = input[4*i+3]
#
#   unscaled_DFT (N/4, input_1, output_1)
#   unscaled_DFT (N/4, input_3, output_3)
#
#   for i in range(N/4):
#	output_1[i] = output_1[i] * W (i, N)
#	output_3[i] = output_3[i] * W (3*i, N)
#
#   for i in range(N/4):
#	tmp_2[i] = output_1[i] + output_3[i]
#	tmp_3[i] = output_1[i] - output_3[i]
#
#   for i in range(N/4):
#	output[i] = tmp_0[i] + tmp_2[i]
#	output[i+N/4] = tmp_1[i] - 1j * tmp_3[i]
#	output[i+N/2] = tmp_0[i] - tmp_2[i]
#	output[i+3*N/4] = tmp_1[i] + 1j * tmp_3[i]

# We didnt do anything here, only reorder the operations. But now, look at the
# first part of this function, up to the calculations of tmp0 and tmp1 : this
# is extremely similar to the radix-2 decimation-in-time algorithm ! or more
# precisely, it IS the radix-2 decimation-in-time algorithm, with size N/2,
# applied on a vector representing the even input coefficients, and giving
# an output vector that is the concatenation of tmp0 and tmp1.
# This is important to notice, because this means we can now choose to
# calculate tmp0 and tmp1 using any DFT algorithm that we want, and we know
# that some of them are more efficient than radix-2...

# This leads us directly to the split-radix decimation-in-time algorithm :

def unscaled_DFT_split_radix_time (N, input, output):
    even_input = vector(N/2)
    input_1 = vector(N/4)
    input_3 = vector(N/4)
    even_output = vector(N/2)
    output_1 = vector(N/4)
    output_3 = vector(N/4)
    tmp_0 = vector(N/4)
    tmp_1 = vector(N/4)

    for i in range(N/2):
	even_input[i] = input[2*i]

    for i in range(N/4):
	input_1[i] = input[4*i+1]
	input_3[i] = input[4*i+3]

    unscaled_DFT (N/2, even_input, even_output)
    unscaled_DFT (N/4, input_1, output_1)
    unscaled_DFT (N/4, input_3, output_3)

    for i in range(N/4):
	output_1[i] = output_1[i] * W (i, N)
	output_3[i] = output_3[i] * W (3*i, N)

    for i in range(N/4):
	tmp_0[i] = output_1[i] + output_3[i]
	tmp_1[i] = output_1[i] - output_3[i]

    for i in range(N/4):
	output[i] = even_output[i] + tmp_0[i]
	output[i+N/4] = even_output[i+N/4] - 1j * tmp_1[i]
	output[i+N/2] = even_output[i] - tmp_0[i]
	output[i+3*N/4] = even_output[i+N/4] + 1j * tmp_1[i]

# This function performs one DFT of size N/2 and two of size N/4, followed by
# N/2 complex multiplies and 3*N/2 complex additions.
# (actually W(0,N) = 1 and W(N/4,N) = -j so we can skip a few of these complex
# multiplies... we will skip 2 for i=0. Also for i=N/8 the W(i,N) and W(3*i,N)
# multiplies can be implemented using only 2 real additions and 2 real
# multiplies)


# We can similarly define the split-radix decimation-in-frequency DFT :

def unscaled_DFT_split_radix_freq (N, input, output):
    even_input = vector(N/2)
    input_1 = vector(N/4)
    input_3 = vector(N/4)
    even_output = vector(N/2)
    output_1 = vector(N/4)
    output_3 = vector(N/4)
    tmp_0 = vector(N/4)
    tmp_1 = vector(N/4)

    for i in range(N/2):
	even_input[i] = input[i] + input[i+N/2]

    for i in range(N/4):
	tmp_0[i] = input[i] - input[i+N/2]
	tmp_1[i] = input[i+N/4] - input[i+3*N/4]

    for i in range(N/4):
	input_1[i] = tmp_0[i] - 1j * tmp_1[i]
	input_3[i] = tmp_0[i] + 1j * tmp_1[i]

    for i in range(N/4):
	input_1[i] = input_1[i] * W (i, N)
	input_3[i] = input_3[i] * W (3*i, N)

    unscaled_DFT (N/2, even_input, even_output)
    unscaled_DFT (N/4, input_1, output_1)
    unscaled_DFT (N/4, input_3, output_3)

    for i in range(N/2):
	output[2*i] = even_output[i]

    for i in range(N/4):
	output[4*i+1] = output_1[i]
	output[4*i+3] = output_3[i]

# The complexity is again the same as for the decimation-in-time varient.


# Now let us now summarize our various algorithms for DFT decomposition :

# radix-2 : DFT(N) -> 2*DFT(N/2) using N/2 multiplies and N additions
# radix-4 : DFT(N) -> 4*DFT(N/2) using 3*N/4 multiplies and 2*N additions
# split-radix : DFT(N) -> DFT(N/2) + 2*DFT(N/4) using N/2 muls and 3*N/2 adds

# (we are always speaking of complex multiplies and complex additions... a
# complex addition is implemented with 2 real additions, and a complex
# multiply is implemented with either 2 adds and 4 muls or 3 adds and 3 muls,
# so we will keep a separate count of these)

# If we want to take into account the special values of W(i,N), we can remove
# a few complex multiplies. Supposing N>=16 we can remove :
# radix-2 : remove 2 complex multiplies, simplify 2 others
# radix-4 : remove 4 complex multiplies, simplify 4 others
# split-radix : remove 2 complex multiplies, simplify 2 others