transforms.py

vlc stand 0.1.99 ist sehr einfach

# This gives the following table for the complexity of a complex DFT :
#	N	real additions	real multiplies	complex multiplies
#       1               0               0               0
#       2               4               0               0
#       4              16               0               0
#       8              52               4               0
#      16             136               8               4
#      32             340              20              16
#      64             808              40              52
#     128            1876              84             144
#     256            4264             168             372
#     512            9556             340             912
#    1024           21160             680            2164
#    2048           46420            1364            5008
#    4096          101032            2728           11380
#    8192          218452            5460           25488
#   16384          469672           10920           56436
#   32768         1004884           21844          123792
#   65536         2140840           43688          269428

# If we chose to implement complex multiplies with 3 real muls + 3 real adds,
# then these results are consistent with the table at the end of the
# www.cmlab.csie.ntu.edu.tw DFT tutorial that I mentionned earlier.


# Now another important case for the DFT is the one where the inputs are
# real numbers instead of complex ones. We have to find ways to optimize for
# this important case.

# If the DFT inputs are real-valued, then the DFT outputs have nice properties
# too : output[0] and output[N/2] will be real numbers, and output[N-i] will
# be the conjugate of output[i] for i in 0...N/2-1

# Likewise, if the DFT inputs are purely imaginary numbers, then the DFT
# outputs will have special properties too : output[0] and output[N/2] will be
# purely imaginary, and output[N-i] will be the opposite of the conjugate of
# output[i] for i in 0...N/2-1

# We can use these properties to calculate two real-valued DFT at once :

def two_real_unscaled_DFT (N, input1, input2, output1, output2):
    input = vector(N)
    output = vector(N)

    for i in range(N):
	input[i] = input1[i] + 1j * input2[i]

    unscaled_DFT (N, input, output)

    output1[0] = output[0].real + 0j
    output2[0] = output[0].imag + 0j

    for i in range(N/2)[1:]:
	output1[i] = 0.5 * (output[i] + conjugate(output[N-i]))
	output2[i] = -0.5j * (output[i] - conjugate(output[N-i]))

	output1[N-i] = conjugate(output1[i])
	output2[N-i] = conjugate(output2[i])

    output1[N/2] = output[N/2].real + 0j
    output2[N/2] = output[N/2].imag + 0j

# This routine does a total of N-2 complex additions and N-2 complex
# multiplies by 0.5

# This routine can also be inverted to calculate the IDFT of two vectors at
# once if we know that the outputs will be real-valued.


# If we have only one real-valued DFT calculation to do, we can still cut this
# calculation in several parts using one of the decimate-in-time methods
# (so that the different parts are still real-valued)

# As with complex DFT calculations, the best method is to use a split radix.
# There are a lot of symetries in the DFT outputs that we can exploit to
# reduce the number of operations...

def real_unscaled_DFT_split_radix_time_1 (N, input, output):
    even_input = vector(N/2)
    even_output = vector(N/2)
    input_1 = vector(N/4)
    output_1 = vector(N/4)
    input_3 = 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)
    # this is again a real DFT !
    # we will only use even_output[i] for i in 0 ... N/4 included. we know that
    # even_output[N/2-i] is the conjugate of even_output[i]... also we know
    # that even_output[0] and even_output[N/4] are purely real.

    unscaled_DFT (N/4, input_1, output_1)
    unscaled_DFT (N/4, input_3, output_3)
    # these are real DFTs too... with symetries in the outputs... once again

    tmp_0[0] = output_1[0] + output_3[0]	# real numbers
    tmp_1[0] = output_1[0] - output_3[0]	# real numbers

    tmp__0 = (output_1[N/8] + output_3[N/8]) * sqrt(0.5)	# real numbers
    tmp__1 = (output_1[N/8] - output_3[N/8]) * sqrt(0.5)	# real numbers
    tmp_0[N/8] = tmp__1 - 1j * tmp__0		# real + 1j * real
    tmp_1[N/8] = tmp__0 - 1j * tmp__1		# real + 1j * real

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

	tmp_0[i] = output_1[i] + output_3[i]
	tmp_1[i] = output_1[i] - output_3[i]

	tmp_0[N/4-i] = -1j * conjugate(tmp_1[i])
	tmp_1[N/4-i] = -1j * conjugate(tmp_0[i])

    output[0] = even_output[0] + tmp_0[0]		# real numbers
    output[N/4] = even_output[N/4] - 1j * tmp_1[0]	# real + 1j * real
    output[N/2] = even_output[0] - tmp_0[0]		# real numbers
    output[3*N/4] = even_output[N/4] + 1j * tmp_1[0]	# real + 1j * real

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

	output[N-i] = conjugate(output[i])
	output[3*N/4-i] = conjugate(output[i+N/4])

# This function performs one real DFT of size N/2 and two real DFT of size
# N/4, followed by 6 real additions, 2 real multiplies, 3*N/4-4 complex
# additions and N/4-2 complex multiplies.


# We can also try to combine the two real DFT of size N/4 into a single complex
# DFT :

def real_unscaled_DFT_split_radix_time_2 (N, input, output):
    even_input = vector(N/2)
    even_output = vector(N/2)
    odd_input = vector(N/4)
    odd_output = 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):
	odd_input[i] = input[4*i+1] + 1j * input[4*i+3]

    unscaled_DFT (N/2, even_input, even_output)
    # this is again a real DFT !
    # we will only use even_output[i] for i in 0 ... N/4 included. we know that
    # even_output[N/2-i] is the conjugate of even_output[i]... also we know
    # that even_output[0] and even_output[N/4] are purely real.

    unscaled_DFT (N/4, odd_input, odd_output)
    # but this one is a complex DFT so no special properties here

    output_1 = odd_output[0].real
    output_3 = odd_output[0].imag
    tmp_0[0] = output_1 + output_3	# real numbers
    tmp_1[0] = output_1 - output_3	# real numbers

    output_1 = odd_output[N/8].real
    output_3 = odd_output[N/8].imag
    tmp__0 = (output_1 + output_3) * sqrt(0.5)	# real numbers
    tmp__1 = (output_1 - output_3) * sqrt(0.5)	# real numbers
    tmp_0[N/8] = tmp__1 - 1j * tmp__0		# real + 1j * real
    tmp_1[N/8] = tmp__0 - 1j * tmp__1		# real + 1j * real

    for i in range(N/8)[1:]:
	output_1 = odd_output[i] + conjugate(odd_output[N/4-i])
	output_3 = odd_output[i] - conjugate(odd_output[N/4-i])

	output_1 = output_1 * 0.5 * W (i, N)
	output_3 = output_3 * -0.5j * W (3*i, N)

	tmp_0[i] = output_1 + output_3
	tmp_1[i] = output_1 - output_3

	tmp_0[N/4-i] = -1j * conjugate(tmp_1[i])
	tmp_1[N/4-i] = -1j * conjugate(tmp_0[i])

    output[0] = even_output[0] + tmp_0[0]		# real numbers
    output[N/4] = even_output[N/4] - 1j * tmp_1[0]	# real + 1j * real
    output[N/2] = even_output[0] - tmp_0[0]		# real numbers
    output[3*N/4] = even_output[N/4] + 1j * tmp_1[0]	# real + 1j * real

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

	output[N-i] = conjugate(output[i])
	output[3*N/4-i] = conjugate(output[i+N/4])

# This function performs one real DFT of size N/2 and one complex DFT of size
# N/4, followed by 6 real additions, 2 real multiplies, N-6 complex additions
# and N/4-2 complex multiplies.

# After comparing the performance, it turns out that for real-valued DFT, the
# version of the algorithm that subdivides the calculation into one real
# DFT of size N/2 and two real DFT of size N/4 is the most efficient one.
# The other version gives exactly the same number of multiplies and a few more
# real additions.


# Now we can also try the decimate-in-frequency method for a real-valued DFT.
# Using the split-radix algorithm, and by taking into account the symetries of
# the outputs :

def real_unscaled_DFT_split_radix_freq (N, input, output):
    even_input = vector(N/2)
    input_1 = vector(N/4)
    even_output = vector(N/2)
    output_1 = 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]

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

    unscaled_DFT (N/2, even_input, even_output)
    # This is still a real-valued DFT

    unscaled_DFT (N/4, input_1, output_1)
    # But that one is a complex-valued DFT

    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[N-1-4*i] = conjugate(output_1[i])

# I think this implementation is much more elegant than the decimate-in-time
# version ! It looks very much like the complex-valued version, all we had to
# do was remove one of the complex-valued internal DFT calls because we could
# deduce the outputs by using the symetries of the problem.

# As for performance, we did N real additions, N/4 complex multiplies (a bit
# less actually, because W(0,N) = 1 and W(N/8,N) is a "simple" multiply), then
# one real DFT of size N/2 and one complex DFT of size N/4.

# It turns out that even if the methods are so different, the number of
# operations is exactly the same as for the best of the two decimation-in-time
# methods that we tried.


# This gives us the following performance for real-valued DFT :
#	N	real additions	real multiplies	complex multiplies
#       2               2               0               0
#       4               6               0               0
#       8              20               2               0
#      16              54               4               2
#      32             140              10               8
#      64             342              20              26
#     128             812              42              72
#     256            1878              84             186
#     512            4268             170             456
#    1024            9558             340            1082
#    2048           21164             682            2504
#    4096           46422            1364            5690
#    8192          101036            2730           12744
#   16384          218454            5460           28218
#   32768          469676           10922           61896
#   65536         1004886           21844          134714


# As an example, this is an implementation of the real-valued DFT8 :

def DFT8 (input, output):
    even_0 = input[0] + input[4]
    even_1 = input[1] + input[5]
    even_2 = input[2] + input[6]
    even_3 = input[3] + input[7]

    tmp_0 = even_0 + even_2
    tmp_1 = even_0 - even_2
    tmp_2 = even_1 + even_3
    tmp_3 = even_1 - even_3

    output[0] = tmp_0 + tmp_2
    output[2] = tmp_1 - 1j * tmp_3
    output[4] = tmp_0 - tmp_2

    odd_0_r = input[0] - input[4]
    odd_0_i = input[2] - input[6]

    tmp_0 = input[1] - input[5]
    tmp_1 = input[3] - input[7]
    odd_1_r = (tmp_0 - tmp_1) * sqrt(0.5)
    odd_1_i = (tmp_0 + tmp_1) * sqrt(0.5)

    output[1] = (odd_0_r + odd_1_r) - 1j * (odd_0_i + odd_1_i)
    output[5] = (odd_0_r - odd_1_r) - 1j * (odd_0_i - odd_1_i)

    output[3] = conjugate(output[5])
    output[6] = conjugate(output[2])
    output[7] = conjugate(output[1])


# Also a basic implementation of the real-valued DFT4 :

def DFT4 (input, output):
    tmp_0 = input[0] + input[2]
    tmp_1 = input[0] - input[2]
    tmp_2 = input[1] + input[3]
    tmp_3 = input[1] - input[3]

    output[0] = tmp_0 + tmp_2
    output[1] = tmp_1 - 1j * tmp_3
    output[2] = tmp_0 - tmp_2
    output[3] = tmp_1 + 1j * tmp_3


# A similar idea might be used to calculate only the real part of the output
# of a complex DFT : we take an DFT algorithm for real inputs and complex
# outputs and we simply reverse it. The resulting algorithm will only work
# with inputs that satisfy the conjugaison rule (input[i] is the conjugate of
# input[N-i]) so we can do a first pass to modify the input so that it follows
# this rule. An example implementation is as follows (adapted from the
# unscaled_DFT_split_radix_time algorithm) :

def complex2real_unscaled_DFT_split_radix_time (N, input, output):
    even_input = vector(N/2)
    input_1 = vector(N/4)
    even_output = vector(N/2)
    output_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] + conjugate(input[N-1-4*i])

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

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

    for i in range(N/4):
	output[i] = even_output[i] + output_1[i].real
	output[i+N/4] = even_output[i+N/4] + output_1[i].imag
	output[i+N/2] = even_output[i] - output_1[i].real
	output[i+3*N/4] = even_output[i+N/4] - output_1[i].imag

# This algorithm does N/4 complex additions, N/4-1 complex multiplies
# (including one "simple" multiply for i=N/8), N real additions, one
# "complex-to-real" DFT of size N/2, and one complex DFT of size N/4.
# Also, in the complex DFT of size N/4, we do not care about the imaginary
# part of output_1[0], which in practice allows us to save one real addition.

# This gives us the following performance for complex DFT with real outputs :
#	N	real additions	real multiplies	complex multiplies
#       1               0               0               0
#       2               2               0               0
#       4               8               0               0
#       8              25               2               0
#      16              66               4               2
#      32             167              10               8
#      64             400              20              26
#     128             933              42              72
#     256            2126              84             186
#     512            4771             170             456
#    1024           10572             340            1082
#    2048           23201             682            2504
#    4096           50506            1364            5690
#    8192          109215            2730           12744
#   16384          234824            5460           28218
#   32768          502429           10922           61896
#   65536         1070406           21844          134714


# Now let's talk about the DCT algorithm. The canonical definition for it is