Faster Mandelbrot image generation using numpy in Python

#
# fractal.py - Mandelbrot image generation using numpy
# Written by Ted Burke - last updated 2-11-2014
#

import numpy

# Define image size and region of the complex plane
w,h = 1280,1280
remin,remax = -2.0, 2.0
immin,immax = -2.0, 2.0

# Create numpy arrays for pixels and complex values
p = numpy.zeros((h, w), dtype=numpy.uint8)
z = numpy.zeros((h, w), dtype=complex)
c = numpy.linspace(remin,remax,w) * numpy.ones((h,w)) + \
    1j * numpy.linspace(immax,immin,h).reshape(h,1) * numpy.ones((h,w))

# Iterative pixel calculation
for n in range(255):
    z = z*z + c
    p = p + (abs(z) < 2.0)

# Output image to binary PGM file
f = open('fractal.pgm', 'wb')
f.write('P5\n{0} {1}\n255\n'.format(w, h))
f.write(p)
f.close()

I ran the program and converted the image to PNG format using ImageMagick as follows:

python fractal.py
convert fractal.pgm fractal.png

fractal

Posted in Uncategorized | Leave a comment

Fractal variations using Python

Every once in a while, I spend an hour or two playing with code to generate fractal images. I find these patterns intriguing and the process of working with them can make for a great workout in complex number arithmetic. When I’m playing with code like this, I find it particularly handy to use Python because it handles the complex numbers very neatly and you can produce a good looking fractal with just a few lines of code. For example, here’s a quick Mandelbrot example:

import numpy

print 'P2\n800 800\n255'

for y in numpy.linspace(-2.0,2.0,800):
    for x in numpy.linspace(-2.0,2.0,800):
        c = x + y*(0+1j)
        z = 0
        for p in range(256):
            z = z*z - c
            if numpy.abs(z) > 4:
                break
        print p,
    print

Please excuse the lack of explanatory comments, but I’m trying to make the point that you don’t need many lines of code!

I ran the program and piped the output into a plain text pgm file, then converted it to PNG format using ImageMagick’s convert command, as shown below:

python mandelbrot.py > mandelbrot.pgm
convert mandelbrot.pgm mandelbrot.png

Ok, here’s the image it produced:

mandelbrot

The process of generating each pixel value in an image like the one above involves starting with a complex number determined by the position within the image (the image spans a rectangular region of the complex plane) and iterating a function recursively to generate a sequence of other points on the complex plane. The faster this iterative process “blows up” (by which I mean that the complex numbers get very large) the darker the pixel.

Today, I was playing with this example and I decided to try to visualise something a little different – for each pixel, I repeated the iteration until the magnitude of the complex value passed a threshold (4 as it happens) then chose the pixel colour based on two things: the distance between the last two complex values and the difference in angle between them. I was thinking of this as “how fast the point is moving” and “how fast the point is revolving around the origin”. Neither of these descriptions are really completely appropriate for an iterative process like this which by definition moves in discrete steps, but that’s how I was thinking of it.

Anyway, the image I produced gave me much food for thought, so I figured I’d post it here in case I want to come back to it again. First, here’s the Python code:

import numpy

print 'P3\n1280 1280\n255'

for y in numpy.linspace(0.0,0.25,1280):
    for x in numpy.linspace(-2.5,-1.25,1280):
        c = x + y*(0+1j)
        z = 0
        for p in range(256):
            newz = z*z + c
            if numpy.abs(newz) > 4:
                z = newz - z
                break
            else:
                z = newz
        r = min(255, int(10 * numpy.abs(z)))
        g = min(255, 128 + int(10 * numpy.angle(z)))
        b = min(255, 128 - int(10 * numpy.angle(z)))
        print r, g, b,
    print

I ran the program and piped the output into a plain text ppm file, then converted it to PNG format using ImageMagick’s convert command, as shown below:

python colours.py > colours.ppm
convert colours.ppm colours.png

This image shows the region of the complex plane where the real part is between -3 and +3 and the imaginary part is between -3j and +3j.

colours1

This image shows the region of the complex plane where the real part is between -1.5 and -1.25 and the imaginary part is between 0j and 0.25j.

colours3

Posted in Uncategorized | Leave a comment

Using SendInput to type unicode characters

I received a query from reader Partha D about generating unicode keystrokes using the SendInput function in Windows. As I understand it, Partha wants to generate one or more unicode keystrokes when a particular keyboard shortcut is pressed. The following example program illustrates the use of the SendInput function to generate keyboard events for unicode characters. It just generates a single keystroke (a Bengali character) after a 3-second delay.

//
// unicodeinput.c - sends a unicode character to whichever
// window is in focus after a delay of 3 seconds (to allow
// time to switch to e.g. Notepad from the command window).
//
// Written by Ted Burke - last updated 2-10-2014
//
// To compile with MinGW:
//
//      gcc unicodeinput.c -o unicodeinput.exe
//
// To run the program:
//
//      unicodeinput.exe
//

// Because the SendInput function is only supported in
// Windows 2000 and later, WINVER needs to be set as
// follows so that SendInput gets defined when windows.h
// is included below.
#define WINVER 0x0500
#include <windows.h>

int main()
{
    // Create a keyboard event structure
    INPUT ip;
    ip.type = INPUT_KEYBOARD;
    ip.ki.time = 0;
    ip.ki.dwExtraInfo = 0;

    // 3 second delay to switch to another window
    Sleep(3000);

    // Press a unicode "key"
    ip.ki.dwFlags = KEYEVENTF_UNICODE;
    ip.ki.wVk = 0;
    ip.ki.wScan = 0x09A4; // This is a Bengali unicode character
    SendInput(1, &ip, sizeof(INPUT));

    // Release key
    ip.ki.dwFlags = KEYEVENTF_UNICODE | KEYEVENTF_KEYUP;
    SendInput(1, &ip, sizeof(INPUT));

    // The End
    return 0;
}

To test the program, I opened a command window and a Notepad window. I compiled and ran the program in the command window, as shown below:

Compiling and running my unicode input example program in a command window

Then I immediately switched the focus to a Notepad window. After three seconds, the following unicode character was typed (I’ve increased the font size to the maximum for clarity):

Screenshot of a Notepad window in which a unicode character has been "typed" by my unicodeinput example program

Posted in Uncategorized | Tagged , , , , , , , , , , , , , , , , | 2 Comments

Three PWM outputs with three different frequencies using the dsPIC30F4011 microcontroller

Arising out of an interesting online discussion with Jaime Mora in the Costa Rica Institute of Technology, I wrote this example program to show how three PWM outputs, each with a different frequency, can be generated using the dsPIC30F4011 microcontroller. Generating three or four PWM outputs with the same frequency is very straightforward, but mutliple PWM outputs with different frequencies is more complicated because each PWM signal must be generated using a different clock source.

Jaime is using the three PWM signals to control IGBT’s in a high-current buck converter. His chosen frequencies are 17 kHz, 19kHz and 23 kHz. My understanding is that these particular frequencies were chosen to avoid harmonic interference. Since 17, 19 and 23 are prime numbers, I think I’m right in saying that the first harmonic frequency that’s common to any two of these PWM signals will be at 323 kHz (because 17 x 19,000 = 19 x 17,000 = 323,000). The three PWM outputs used in this program are:

  • A: PWM1H on pin 37. This 17 kHz PWM signal is generated by the “Motor Control PWM” module. The clock source is the dedicated PWM clock in that module.
  • B: OC1 on pin 23. This 19 kHz PWM signal is generated using Output Compare channel 1 (in PWM mode), with Timer 2 as its clock source.
  • C: OC2 on pin 18. This 23 kHz PWM signal is generated using Output Compare channel 2 (in PWM mode), with Timer 3 as its clock source.

I included a handy function called set_duty_cycles, which allows all three channels’ duty cycles to be set at once. Each duty cycle should be specified as a floating point value between 0 and 1. This is the function prototype:

void set_duty_cycles(float a, float b, float c);
  • a is the duty cycle (between 0 and 1) for PWM1H on pin 37.
  • b is the duty cycle (between 0 and 1) for OC1 on pin 23.
  • c is the duty cycle (between 0 and 1) for OC2 on pin 18.

I’ve tested this example program with LEDs attached to the three output pins so that I could confirm that PWM was working and that the duty cycle really was varying on each channel (as shown in the video below). Jaime was kind enough to send me the following screen captures from his oscilloscope when he was testing his circuit. They seem to confirm that the PWM frequencies are 17 kHz, 19 kHz and 23 kHz, as expected.

Here’s a video of the three LEDs pulsating:

Here’s the complete source code:

//
// This example shows how the dsPIC30F4011 can generate three PWM
// signals, each with a different frequency (17000, 19000, 23000 Hz)
//
// To generate three different frequencies, three separate clock sources
// must be used, which necessitates some mixing and matching:
//
//   PWM signal A (17 kHz) is generated using the PWM module with its
//   own timebase (PTMR) used as the clock source.
//
//   PWM signal B (19 kHz) is generated using Output Compare channel 1
//   with Timer 2 as the clock source.
//
//   PWM signal C (23 kHz) is generated using Output Compare channel 2
//   with Timer 2 as the clock source.
//
// I compiled this with Microchip's XC16 C compiler, using the
// following commands:
//
//   xc16-gcc main.c -mcpu=30F4011 -Wl,--script=p30F4011.gld
//   xc16-bin2hex a.out
//
// Written by Ted Burke, Last updated 24-9-2014
//
 
#include <xc.h>
#include <libpic30.h>
#include <math.h>
 
// Configuration settings
_FOSC(CSW_FSCM_OFF & FRC_PLL16); // Fosc=16x7.5MHz, i.e. 30 MIPS
_FWDT(WDT_OFF);                  // Watchdog timer off
_FBORPOR(MCLR_DIS);              // Disable reset pin
 
// Function prototype
void set_duty_cycles(float a, float b, float c);
 
int main(void)
{
    // Use PWM module for PWM output A with PTMR as clock source
    PWMCON1 = 0x0011;       // Enable PWM1 (high and low pins)
    PTCONbits.PTCKPS = 0;   // prescale=1:1 (0=1:1, 1=1:4, 2=1:16, 3=1:64)
    PTPER = 1764;           // 17 kHz PWM frequency (PTPER + 1 = 30000000 / 17000)
    PDC1 = 1200;            // 50% duty cycle on PWM channel 1
    PTMR = 0;               // Clear 15-bit PWM timer counter
    PTCONbits.PTEN = 1;     // Enable PWM time base
 
    // Use OC1 for PWM output B with Timer 2 as clock source
    PR2 = 1578;             // 19 kHz PWM frequency (PR2 + 1 = 30000000 / 19000)
    OC1R = PR2 / 2;         //
    OC1RS = PR2 / 2;        // Select 50% duty cycle initially
    OC1CONbits.OCTSEL = 0;  // Select Timer 2 as clock source for OC1
    OC1CONbits.OCM = 0b110; // PWM mode
    T2CONbits.TON = 1;      // Turn on Timer 2
 
    // Use OC2 for PWM output C with Timer 3 as clock source
    PR3 = 1303;             // 23 kHz PWM frequency (PR3 + 1 = 30000000 / 23000)
    OC2R = PR3 / 2;         //
    OC2RS = PR3 / 2;        // Select 50% duty cycle initially
    OC2CONbits.OCTSEL = 1;  // Select Timer 3 as clock source for OC2
    OC2CONbits.OCM = 0b110; // PWM mode
    T3CONbits.TON = 1;      // Turn on Timer 3
    
    // Make OC1 and OC2 outputs (same pins as RD0 and RD1)
    TRISD = 0b1111111111111100;
    
    //
    // Now, vary duty cycles on the three PWM outputs so that the
    // three LEDs pulsate.
    //
    // I'm using an inverted squared sine waveform to pulsate
    // each LED. Between each pair of LEDs, there's a phase shift
    // of 2*pi/3 radians.
    //
    // It doesn't really matter what's actually happening in this part
    // of the program - it's just something to show that the duty
    // cycle really is variable on each PWM output.
    //
    float t=0, pi=3.1428, f=0.25, a, b, c;
    while(1)
    {
        a = sin(2 * pi * f * t);
        a = 1 - a * a;
        b = sin(2 * pi * f * t + 2 * pi / 3);
        b = 1 - b * b;
        c = sin(2 * pi * f * t + 4 * pi / 3);
        c = 1 - c * c;
        
        set_duty_cycles(a, b, c);
        
        __delay32(300000); // 10 ms delay
        t += 0.01; // t is in seconds
    }
}
 
//
// This function provides a way to set the duty cycle on all three
// PWM output channels in one go. Arguments a, b and c are the three
// desired duty cycle values - each should be in the range from 0 to 1.
//
void set_duty_cycles(float a, float b, float c)
{
    PDC1 = a * 2 * PTPER;
    OC1RS = b * PR2;
    OC2RS = c * PR3;
}
Posted in Uncategorized | Leave a comment

8-channel PWM with the MSP430G2553

Last week, I spent a bit of time working out how to generate PWM for servo control using the MSP430G2553. It took me a couple of hours to work it out, but worked well once I got it going. Today, I’ve been experimenting with a small robot manipulator that we’re exhibiting in our open day tomorrow in DIT Kevin St. This manipulator has a total of seven servo motors, but only 6 degrees of freedom because two of the servos are doubled up on the shoulder joint and are controlled by the same PWM signal. I wanted to control this manipulator with an MSP430G2553, but it doesn’t provide that many dedicated PWM outputs, so I needed to investigate an alternative way of doing it. The code below is my ad hoc solution.

//
// 8-channel PWM example for MSP430G2553
//
// Written by Ted Burke - last updated 4-4-2014
//

#include <msp430.h>

// PWM channels duty cycle array
int pw[] = {1000,1000,1000,1000,1000,1000,1000,1000};

int main( void )
{
    WDTCTL = WDTPW + WDTHOLD; // Disable watchdog timer

    P1DIR = 0b00111111; // Make P1.0-5 outputs

    // Configure Timer A0 Compare interrupts
    TA0CTL = TASSEL_2 + MC_1 + ID_0; // "up" mode, input divider = 1
    TA0CCR0 = 2500;                  // set timer period to 2.5ms
    TA0CCTL0 = CCIE;                 // enable CC0 interrupt
    TA0CCR1 = 1500;                  // set pulse width to 1.5ms
    TA0CCTL1 = CCIE;                 // enable CC1 interrupt
    _BIS_SR(GIE);                    // global interrupt enable

    // From this point on, we only need to write values
    // into the pw[] array to set the duty cycle on all
    // eight PWM channels (P1.0-7), or to be precise,
    // whichever channels are actually enabled as digital
    // outputs (six in this case).

    //
    // A quick example to test: Do a different number of
    // angle steps on each of the six PWM outputs. 1 step
    // on channel 0, 2 steps on channel 1, 3 steps on
    // channel 2, and so on.
    //
    int channel, counter = 0;
    while(1)
    {
        counter++;
        for (channel=0 ; channel<6 ; ++channel)
        {
            pw[channel] = 1000 + (counter%(channel+1))*100;
        }
        __delay_cycles(500000);
    }

    return 0;
}

//
// Timer A0 CC0 interrupt service routine.
// This ISR is triggered when Timer A0 reaches
// its maximum value TA0CCR0 and resets to zero.
// This ISR starts a pulse on one PWM channel
// every 2.5ms. It cycles through all eight channels
// in turn. After starting the pulse, it sets
// the TA0CCR1 register to the pulse width of
// the current pin so that the pulse will be
// ended by the partner ISR after the appropriate
// time delay (for this channel's duty cycle).
//
#pragma vector=TIMER0_A0_VECTOR
__interrupt void Timer_A0_CC0(void)
{
    static n=0;         // PWM channel index

    P1OUT = 1 << n;     // Set P1.n high
    TA0CCR1 = pw[n];    // Set timer for current pin's pulse width

    n = (n+1)%8;        // Move on to next PWM channel
    TA0CCTL0 &= ~CCIFG; // Reset CC interrupt flag
}

//
// Timer A0 CC1 interrupt service routine.
// This ISR is responsible for ending each PWM
// pulse on all channels. It is triggered when
// Timer A0 matches TA0CCR1, which is at a
// different time for each PWM channel.
//
#pragma vector=TIMER0_A1_VECTOR
__interrupt void Timer_A1_CC1(void)
{
    P1OUT = 0;
    TA0CCTL1 &= ~CCIFG;
}

The following video shows what the above program does. In the video you can see me switching the signal controlling the servo motor from one PWM output to the next in the following sequence: P1.0, P1.1, P1.2, P1.3, P1.4, P1.5.

This method involves cooperation between two Timer A0 interrupt service routines (ISRs), one which is responsible for starting pulses on every channel (TA0 CC0) and one which is responsible for ending pulses on every channels (TA1 CC1). The pulses for the six channels are interleaved, in the sense that the channels take it in turns to output one pulse. There is never more than one PWM channel outputting a pulse at any time. Even at its longest pulse width, the servo control PWM has a relatively low duty cycle (approximately 10%), so there is no difficulty interleaving 8 channels in this fashion, as illustrated below:

8-channel_pwm

Download Editable SVG version of above image

Posted in Uncategorized | 2 Comments

Simple Phaser Framework example – Flappy Words

I spent some time this week playing Aoife Crowley’s Flappy Bod game (my top score = 8). I was really impressed by the simplicity of the Phaser Framework which it runs on, so I decided to try it out for myself. Phaser is an open source framework for 2-D web game development. It’s really easy to get started with and I was able to get a few sprites moving around on the screen and responding to user input within minutes.

Of course, I still have to decide exactly how my own Flappy Words game is going to work, but the basic idea is to facilitate text communication (i.e. spelling messages out) using a single switch (the space bar), in something like the style of Flappy Bird. I suppose I’m picturing a cross between Flappy Bird and Dasher. The exact details remain to be worked out, but I thought I might as well document my program as it currently stands, which is a very basic indeed – no spelling at all yet; just one flying teapot.

Click on the screenshot below to see what this example does so far (not a lot!):
Flappy Words Screenshot - 03222014 - 01:53:13 AM

The game basically consists of just a few files:

  1. index.html – This is the main HTML file that you actually point the browser at. It doesn’t contain much apart from a div with the same id (“flappy-words”) as that specified when the game object is created in the game.js JavaScript File. The contents of this file are shown below.
  2. game.js – This is the file that contains all my JavaScript code for the game so far. The contents of this file are shown below.
  3. background.png – This is the background wallpaper used in the game world.
  4. teapot.png – This is the PNG image used for the flapping teapot sprite.
  5. phaser.min.js – This file contains the entire Phaser Framework. This one is version 2.0 and it’s exactly as I downloaded it from the Phaser github repo.

While developing with Phaser, of course you need to be able to test the program repeatedly in a browser. For security reasons, browsers will not generally allow javascript code to access the local filesystem, so to test the game while you’re writing it you need to set up a web server and make the web browser access the files through it. The easiest way to do this is to use the built in web server in the Python interpreter. I only discovered this brilliant feature yesterday in an article from the Linux Journal. Is there nothing Python can’t do?!

I store all my game files in the same folder and then start Python in that folder using the following command:

python -m SimpleHTTPServer

By default, the Python web server accepts connections on port 8000 (although you can specify a different port if necessary). I therefore test my game by pointing Firefox at the following URL:


http://127.0.0.1:8000

By the way, this is the Inkscape SVG file that I created to make the PNGs for the teapot sprite and background:

  1. teapot.svg

This is the HTML code in “index.html”:

<!DOCTYPE html>

<!-- This is the main HTML file for Flappy Words -->
<!-- Written by Ted Burke, last updated 21-3-2014 -->

<html>

<head>
	<title>Flappy Words</title>
	
	<style>
		#flappy-words {
			width: 800px;
			margin: auto;
			margin-top: 10px;
		}
	</style>
	
	<script type="text/javascript" src="phaser.min.js"></script>
	<script type="text/javascript" src="game.js"></script>
	
</head>

<body>
	<h1 style="text-align: center;">Flappy Words</h1>
	<div id="flappy-words"></div>
</body>

</html>

This is the JavaScript code in “game.js”:

//
// game.js - Main JavaScript file for Flappy Words
// Written by Ted Burke - last updated 21-3-2014
//

var game = new Phaser.Game(800, 600, Phaser.AUTO, 'flappy-words', {preload: preload, create: create, update: update});

var teapot;
var text;
var counter = 0;

function preload () {
	game.load.image('teapot', 'teapot.png');
	game.load.image('background', 'background.png');
}

function create() {
	game.add.image(0, 0, 'background');
	
	teapot = game.add.sprite(game.world.centerX, game.world.centerY, 'teapot');
	teapot.inputEnabled = true;
	teapot.events.onInputDown.add(flap, this);
	teapot.scale.x = 0.4;
	teapot.scale.y = 0.4;
	game.physics.enable(teapot, Phaser.Physics.ARCADE);
	teapot.anchor.setTo(0.5,0.5);
	teapot.body.collideWorldBounds = true;
	teapot.body.gravity.y = 200;
	teapot.body.bounce.y = 0.4;
	teapot.body.maxAngular = 200;
	//teapot.body.angularDrag = 20;
	
	text = game.add.text(250, 16, 'Press SPACE to flap...', { fill: '#ffffff' });
	
	spaceKey = game.input.keyboard.addKey(Phaser.Keyboard.SPACEBAR);
	spaceKey.onDown.add(flap, this);
}

function update() {
	if (spaceKey.isDown)
	{
		text.fill = '#000000';
		teapot.body.angularVelocity = 100;
	}
	else
	{
		text.fill = '#ffffff';
		teapot.body.angularAcceleration = -50.0 * teapot.angle - 5.0 * teapot.body.angularVelocity;
	}
}

function flap() {
	teapot.body.velocity.y = -200;
	counter++;
	text.text = "You flapped " + counter + " times!";
}

These are the sprite and background PNG images:

Posted in Uncategorized | Leave a comment

Installing Microchip XC16 in CrunchBang Linux

XC16 is Microchip’s C compiler for the 16-bit dsPIC microcontrollers (PIC24, dsPIC30F, dsPIC33). I sometimes need to use it when I’m working in a live session in CrunchBang Linux.

First download the Linux installer for XC16

Then change directory into the downloads folder (/home/crunchbang/downloads/ in my case) and extract the installer from the tar archive:

cd downloads
tar xvf xc16-v1.21-linux-installer.run.tar

Now, log in as root, install some dependencies, and run the installer:

sudo su
apt-get update
apt-get install lib32z1
./xc16-v1.21-linux-installer.run

You will be asked to agree to a software license and to answer a few questions on installation options. The default options are mostly acceptable, but make sure you reject the pay options and time-limited premium license. The free time-unlimited version should be adequate for most purposes.

By default the compiler is installed in /opt/microchip/xc16/

To compile a program, here’s an example command line:

/opt/microchip/xc16/v1.21/bin/xc16-gcc main.c -mcpu=30F4011 -Wl,--script=p30F4011.gld

That produces a file called a.out. To convert the compiled program to hex format (i.e. “a.hex”), use the bin2hex tool, as follows:

/opt/microchip/xc16/v1.21/bin/xc16-bin2hex a.out
Posted in Uncategorized | Leave a comment