Complex Numbers

[padx]

Hi everyone

I tried to implement a fft algorithm on my pyboard. After a few workarounds i actually succeeded. But I soon found the problem that working with complex numbers doesn't seem to be completed. For example I cannot get the real or imaginary part. So my question is if someone is currently working on complex functions or how I can do it by my own.
I tried to read through the different tutorials but I did not yet found one for extending existing types.

Hope anyone can help.

[blmorris]

You are right that the .real and .imag attributes of complex numbers aren't implemented yet; I also noticed this while trying to implement FFT; unfortunately I didn't pursue it further because my project had more pressing priorities at the time. I assume that it it a straightforward addition to 'py/objcomplex.c' but my minimal skills in C aren't up to the task. Probably should have posted it as an issue on GitHub: https://github.com/micropython/micropyt ... state=open but I let it slide instead

I'm curious how you managed to implement FFT in spite of the missing attributes - would you mind posting your code?

-Bryan

[blmorris]

Possibly of interest- another user on the forum started documenting his experience implementing a module in C: http://forum.micropython.org/viewtopic.php?f=3&t=176
He is attempting to implement the CAN bus on the pyboard; maybe not directly relevant to this, but has some explanations on the micro python code structure that I haven't seen elsewhere.

-Bryan

[padx]

Well actually for the FFT itself you do not need the imaginary or real part at least as far as I understand it^^
so here is my code (is not the fastest but was my first try without numpy or similiar stuff and actually I don't know if it is actually working as it should. I get a lot of noise):
Here the fft function:

Code: Select all

def fft(x):
	n = len(x)
	if n == 1:
		return x
	else:
		xeven = [x[i] for i in range(0, n, 2)]
		xodd = [x[i] for i in range(1,n,2)]
		#print(xeven)
		Feven = fft(xeven)
		Fodd = fft(xodd)
		#print(Feven)
		combined = [0] *n
		for m in range(0,int(n/2)):
			c = math.e**((-2.0j *math.pi*m)/n) *Fodd[m]
			combined[m] = Feven[m] +  c
			combined[m + int((n-1)/2)]= Feven[m] - c
		return combined

And here my complete code:

Code: Select all

import pyb
import math
import cmath

lcd = pyb.LCD('Y')
lcd.light(True)

def printValuesOnScreen(values, offset=0):
	x = 0
	xcoord = offset
	while xcoord < 200:
		#print(x)
		if x < len(values):
			#print(values[x]/10)
			xval = 0
			for i in range(0,1):
				xval += int(values[x+i]/10)
			lcd.pixel(xcoord,27 - int(xval) ,1)
			#lcd.pixel(xcoord,54 - int(xval),1)
		lcd.pixel(xcoord,15,1)
		x+=1
		xcoord += 1
	lcd.show()

def printDFFTOnScreen(values, width=30):
	x = 0
	xcoord = 0
	while xcoord < width:
		#print(x)
		if x < len(values):
			#print(values[x]/10)
			xval = 0
			for i in range(0,1):
				xval += int(values[x+i]/5)
			#lcd.pixel(xcoord,27 - int(xval) ,1)
			for y in range(int(xval)):
				lcd.pixel(xcoord,30 - y,1)
		lcd.pixel(xcoord,30,1)
		x+=1
		xcoord += 1
	lcd.show()

def fft(x):
	n = len(x)
	if n == 1:
		return [x[0]/10]
	else:
		xeven = [x[i] for i in range(0, n, 2)]
		xodd = [x[i] for i in range(1,n,2)]
		#print(xeven)
		Feven = fft(xeven)
		Fodd = fft(xodd)
		#print(Feven)
		combined = [0] *n
		for m in range(0,int(n/2)):
			c = math.e**((-2.0j *math.pi*m)/n) *Fodd[m]
			combined[m] = Feven[m] +  c
			combined[m + int((n-1)/2)]= Feven[m] - c
		return combined
	

adc = pyb.ADC('X22')
while(True):
	buf = bytearray(200)
	adc.read_timed(buf, 640000)
	lcd.fill(0)
	coefs = fft([buf[i] for i in range(0,180,6)])
	arr = [abs(coefs[i]) for i in range(len(coefs)-1)]
	#print(arr)
	printDFFTOnScreen(arr)
	printValuesOnScreen(buf, 30)
	

I'm playing around with the microphone and trying to use the fft for frequency analysis. I also had to use the list comprehension because sclices with stepsizes more than one are also not implemented.
But I could not yet flash the newest rom, so I'm waiting with any bug report...

Patrick

[pfalcon]

I also had to use the list comprehension because sclices with stepsizes more than one are also not implemented.
But I could not yet flash the newest rom

Yes:

Code: Select all

$ micropython 
Micro Python v1.1.1-141-g2cf3810 on 2014-07-12; linux version
>>> [1,2,3,4][::2]
[1, 3]
>>> 

[padx]

Ok so I'll definitely have to flash a newer rom

[pythoncoder]

At the risk of stating the obvious you should be able to get the real and imaginary parts with a little maths. Bear in mind I haven't yet got a MicroPython board, so it's possible that I'm using un-implemented functionality here, but this works in standard Python.

Code: Select all

import cmath
from math import cos,sin
z = 3 +7j
r, phi = cmath.polar(z)
re, imag = r*cost(phi), r*sin(phi)

Regards, Pete

[pythoncoder]

I was struck by the simplicity of this code compared to a routine I've used. Apologies if I'm misunderstanding its application, but my tests seem to produce some odd results. I used the following to test the fft routine:

Code: Select all

import math

def test(length, sine =False):
   if sine:
      a = [math.sin(2*math.pi*n/length) for n in range(length)]
      print("Sine samples")
      for z in a:
         print(z)
      b = fft(a)
      print("Sine Results")
      for z in b:
         print(abs(z))
   else:
      a  = [1.0]*(length//2) + [-1.0]*(length//2)
      if length % 2:
         a += [-1.0] # cope with odd number of samples
      print("Square samples")
      for z in a:
         print(z)
      b = fft(a)
      print("Square Results")
      for z in b:
         print(abs(z))

def fft(x):
   n = len(x)
   if n == 1:
      return [x[0]/10]
   else:
      xeven = [x[i] for i in range(0, n, 2)]
      xodd = [x[i] for i in range(1,n,2)]
      #print(xeven)
      Feven = fft(xeven)
      Fodd = fft(xodd)
      #print(Feven)
      combined = [0] *n
      for m in range(0,int(n/2)):
         c = math.e**((-2.0j *math.pi*m)/n) *Fodd[m]
         combined[m] = Feven[m] +  c
         combined[m + int((n-1)/2)]= Feven[m] - c
      return combined

Firstly the results returned depend strongly on the length of the list. Results where the length is a power of two are revealing: the square wave results show only a single nonzero frequency bin. Run the test with a length of 63 with sine true and numerous bins become populated. It would seem that something is wrong here, but I can't claim to understand the code to fix it!

Regards, Pete

[padx]

Yeah I've got also odd results. And yes this algorithm is only meant to be used with powers of two, because you are always dividing the array by two, hence if it is not a power of two, the results are strange (btw. thank you for pointing that out, didn't think about that).

I'll post the solution as soon I got time to flash new rom and test everything if I find time I'll also try to implement the missing functions. At the moment I should learn for my exams in august, so it is kind of a bad idea doing too much with the board

Greets

[pythoncoder]

I've investigated this a further and found the following code, which is similar to yours but produces results which look more convincing. As you point out, the array length should be a power of two. The output list contains the frequency components followed by their complex conjugates which can be ignored. Element zero is the DC component which would normally be zero in an audio application. Note that the results require scaling in proportion to the array length, as in my printresults() test function.
Here is the code (tested on a PC):

Code: Select all

def fft(x):
    n = len(x)
    if n <= 1:
        return x
    even = fft(x[0::2])
    odd =  fft(x[1::2])
    return [even[m] + math.e**(-2j*math.pi*m/n)*odd[m] for m in range(n//2)] + \
           [even[m] - math.e**(-2j*math.pi*m/n)*odd[m] for m in range(n//2)]

And here is the code I used to test it

Code: Select all

import math

fformat = "{:10.4f}"

def test(length, sine =False):
   assert checklength(length),"Length of array must be a power of two"
   if sine:
      a = [math.sin(2*math.pi*n/length) for n in range(length)]
      print("Sine samples")
      for z in a:
         print(fformat.format(z))
      b = fft(a)
      printresults(b,"Sine Results")
   else:
      a  = [1.0]*(length//2) + [-1.0]*(length//2)
      print("Square samples")
      for z in a:
         print(z)
      b = fft(a)
      printresults(b, "Square Results")

def printresults(b, msg):   # Print magnitude only. Ignore complex conjugates.
    print(msg)              # scale by array length
    for z in b[0:len(b)//2]:
        print(fformat.format(2*abs(z)/len(b)))  

def checklength(n):         # Sanity check before actually doing the FFT
    while n and not n & 1:
        n >>=1
    return n == 1

# Enter with an array of complex numbers of length 2**N
def fft(x):
    n = len(x)
    if n <= 1:
        return x
    even = fft(x[0::2])
    odd =  fft(x[1::2])
    return [even[m] + math.e**(-2j*math.pi*m/n)*odd[m] for m in range(n//2)] + \
           [even[m] - math.e**(-2j*math.pi*m/n)*odd[m] for m in range(n//2)]

A final comment. While this method is extremely compact, the more complex code I'm familiar with exists for a reason. The compact code is recursive and creates and destroys lists on the stack which is costly in terms of memory use and processor cycles. The conventional algorithms perform the DFT in place, with the results appearing in the array which once contained the samples. They use iteration rather than recursion. Consequently they are more efficient.

In your application with small sample sets this is probably irrelevant.

Please accept my apologies if you know all this already!

Regards, Pete