2

I am trying to implement recursive FFT.

import numpy as np
from math import e, pi

def rdft(a):
    n = a.size
    if n == 1:
        return a
    i = complex(0, 1)
    w_n = e ** (2 * i * pi / float(n))
    w = 1
    a_0 = np.zeros(int(math.ceil(n / 2.0)))
    a_1 = np.zeros(n / 2)
    for index in range(0, n):
        if index % 2 == 0:
            a_0[index / 2] = a[index]
        else:
            a_1[index / 2] = a[index]
    y_0 = rdft(a_0)
    y_1 = rdft(a_1)
    y = np.zeros(n)
    for k in range(0, n / 2):
        y[k] = y_0[k] + w * y_1[k]
        y[k + n / 2] = y_0[k] - w * y_1[k]
        w = w * w_n
return y



if __name__ == "__main__":
    a = np.array([1, 0,0, -2])
    print rdft(a)

It gives me the following error:

[-1.  1.  3.  1.]
/path/file.py:22: ComplexWarning: Casting complex values to real discards the imaginary part y[k] = complex(y_0[k] + w * y_1[k])
/path/file.py:23: ComplexWarning: Casting complex values to real discards the imaginary part y[k + n / 2] = y_0[k] - w * y_1[k]

I'm not familiar with complex calculations in Python. The real values are correct but as it says it discards the imaginary part.

mpourreza
  • 177
  • 1
  • 1
  • 15
  • Maybe you want to explain, what your function is supposed to do. When you say polynomial multiplication, you don't mean [this](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.polymul.html#numpy.polynomial.polynomial.polymul), I assume? And please edit your question and post the full traceback. – Mr. T Jan 31 '18 at 23:19
  • You can write the complex part of a number with a `j` suffix, like `2j * pi` as you use it. – Chiel Jan 31 '18 at 23:29
  • This code does not run. Have you tested yourself what you have pasted? – Chiel Jan 31 '18 at 23:31
  • Yes, it runs and prints the output. But without the imaginary part. – mpourreza Jan 31 '18 at 23:32
  • The first error message doesn't match your code. – Mr. T Jan 31 '18 at 23:33
  • Yes. it adds the complex() in error code but my code doesn't have it. – mpourreza Jan 31 '18 at 23:34
  • Strange. Not in my error message. – Mr. T Jan 31 '18 at 23:44

1 Answers1

2

You try to cast complex numbers to a numpy array that has been defined as float. To overcome this problem, define y as a numpy array with complex values:

y = np.zeros(n, dtype = complex)

Output after this change:

[-1.+0.j  1.+2.j  3.+0.j  1.-2.j]

If you want to know more about numpy data types, read the numpy user guide.

Mr. T
  • 11,960
  • 10
  • 32
  • 54