C or Python? Comparison of execution time for Mandelbrot image generation

Because I’ve been generating a lot of fractal images recently, I’ve found myself spending quite a bit of time sitting around waiting for programs I’ve written to finish long repetitive calculations. In particular, I’ve generated several long sequences of several hundred images that I subsequently combined into animations.

One of the many great things about teaching electrical engineering is that I get to dip in and out of different programming languages all the time, but C remains my primary language and it’s still what I’m most likely to reach for when I’m writing something for my own curiosity. However, over the last couple of months I’ve been using Python for most of my fractal programming because I find its complex number syntax elegant. Furthermore, although complex numbers have been natively supported in C since the C99 standard, I’ve just never really embraced that language feature (until now).

Anyway, because I’ve been spending all this time waiting while numbers are getting crunched, I decided to do a quick comparison of execution time for the same task between Python and C. Naturally, since C is compiled and Python is interpreted, it comes as no surprise that C is significantly faster. Nevertheless, I was still surprised by how much faster it is.

The task I chose for this experiment was the generation of a Mandelbrot Set image of resolution 4000×4000 pixels, spanning a square region of the complex plane from -2 to +2 on the real axis and from -2j to +2j on the imaginary axis. This is the image produced:

mandelbrot

I implemented the same algorithm in C and Python and then performed a rough measurement of the execution time of each program by running the “date” command line tool immediately before and after it, as shown below…

Screenshot_2016-02-02_23-56-08

As you can see, the approximate execution times were as follows:

  • Python implementation: 105 seconds
  • C implementation: 4 seconds

In other words, the C implementation executes more than 26 times faster than the Python implementation! Compiling the program with gcc (also shown in the above screen shot) only takes a fraction of a second too. I’m sure a better programmer could make either of these implementations run a bit faster, but unless the Python version can run dramatically faster than this then the difference is more than enough to persuade me to embrace complex maths in C.

This is the full C code:

//
// mandelbrot.c - written by Ted Burke - 2-Feb-2016
//
// To compile:
//      gcc mandelbrot.c -o mandelbrot -lm
//

#include <complex.h>
#include <stdio.h>

// Create an array of bytes to store pixel values
unsigned char p[4000][4000];

void main()
{
    double complex c, z;
    int w=4000, h=4000, x, y;
    unsigned char px;
    
    // Open output file
    FILE* f = fopen("mandelbrot_c.pgm","w");
    fprintf(f, "P5\n%d %d\n255\n", w, h);
    
    // Calculate pixel values
    for (y=0 ; y<h ; ++y)
    {
        for (x=0 ; x<w ; ++x)
        {
            // Initialise complex values c and z for next
            // iterative pixel value calculation
            c = (-2.0 + x*4.0/w) + (2.0 - y*4.0/h)*I;
            z = 0;
            
            // Iterate complex function while darkening pixel
            px = 255;
            while(px >= 10 && cabs(z) < 4.0)
            {
                z = z*z + c;
                px -= 10;
            }
            
            // Store latest pixel value in byte array
            p[y][x] = px;
        }
    }
    
    // Write the byte array to the output file
    fwrite(p, 1, w*h, f);
    fclose(f);
}

This is the full Python code:

#
# mandelbrot.py - written by Ted Burke - 2-Feb-2016
#

import numpy

# Create an array of bytes to store pixel values
w,h = 4000,4000
p = numpy.zeros((h,w),dtype=numpy.uint8)

# Open output file
f = file('mandelbrot_python.pgm','w')
f.write('P5\n{:d} {:d}\n255\n'.format(w, h))

# Calculate pixel values
for y in range(h):
    for x in range(w):
        # Initialise complex values c and z for next
        # iterative pixel value calculation
        c = (-2.0 + x*4.0/w) + (2.0 - y*4.0/h)*(0+1j)
        z = 0 + 0j
        # Iterate complex function while darkening pixel
        px = 255
        while px >= 10 and abs(z) < 4.0:
            z = z*z + c
            px -= 10
        # Store latest pixel value in byte array
        p[y][x] = px

# Write the byte array to the output file
f.write(bytearray(p))
f.close()
Posted in Uncategorized | 4 Comments

Julia Set Explorer

In the past, writing simple programs to generate fractal images is something I have found myself returning to over and over again. I’ve done it so many times that I can probably write a Mandelbrot Set generator in my sleep. However, I’ve really struggled to understand at an intuitive level why these shapes emerge the way they do from such deceptively simple looking iterative functions.

Over the recent Christmas break I spent some time playing with Julia Set visualisations. My aim was to gain a better understanding of why these sets often take on exotic fractal shapes. I’ve really stuck with it this time and I finally feel like I’m beginning to gain a better understanding. Today, I finally feel like I’ve achieved a bit of a breakthrough in my own understanding, which I’ll try to illustrate in the following images.

The solid black shape in each image is the so-called filled Julia Set of a complex iterating function, f. Specifically, in the images above,

f(z) = z^2 + c

where z and c are complex values and c = 0.325 + 0.285j. Different values of c produce very different shapes. The basic idea is that for a given value of c you select an initial value of z – let’s call it z0 – and then apply the iterating function f over and over again to generate a sequence of complex values

z_0, z_1, z_2, z_3,...

where

z_{n+1}=f(z_n) \textrm{  for  } n>0

We call this sequence the orbit of z0 under f. For the function f shown above, every point in the complex plane falls into one of three possible categories:

  • The points which are coloured white in the images above have orbits which spiral out to infinity. These points make up F0, one of the Fatou domains of f.
  • The points which are coloured black in the images above have orbits which remain bounded (they stay within the black region). These points make up F1, the second Fatou domain of f.
  • The third category is the set of points which form the boundary between the two Fatou domains – this is the Julia Set of f. In the images above, the Julia Set is basically where the black region meets the white region. The orbit of any point in the Julia Set lies completely within the Julia Set – i.e. the orbit just jumps around the boundary forever, never going inside or outside.

The Julia Set of this function has a fractal shape, as is the case for many complex functions. What I’m really trying to understand – in an intuitive sense – is why that is the case. Self-similarity is a frequently mentioned characteristic of fractals – specifically the recurrence of similar shapes within the fractal at different levels of magnification. There’s no doubt that very similarly shaped outcrops appear repeatedly along the boundary of the filled Julia sets shown above. Furthermore, if we were to zoom into that boundary, we would see identically shaped outcrops recurring over and over again at different levels of magnification.

The insight that struck me so forcefully today relates to picturing how the iterating function generates this self-similarity. In the image on the left, the green dot marks the location of z0 in the complex plane. The orbit of z0 is shown as a series of red dots joined by blue lines which illustrate what happens to z each time the function is iterated (first z is squared which is represented by the curve, then c is added to it which is represented by the straight line). When I traced z0 along just outside the boundary, I could see z1, z2, z3, etc tracing out similar shapes further up the boundary. The best way I can describe it is that it reminded me of a pantograph, but with many pens rather than just one.

The images above are screen shots of a python program (“julia.py”) I wrote to display filled Julia sets and visualise the orbits of different complex numbers. The complete Python code is shown below. To run this, I think you just need Python and Tkinter. I’m running Python 2.7.9, but hopefully it should work in Python 3 also (?). Basic instructions are below the code listing.

#
# julia.py - Real time orbit viewer for iterative complex functions
# written by Ted Burke, last updated 21-1-2016
#

from Tkinter import *
import math
import sys

def generate_fractal():
    sys.stdout.write('Generating fractal for c = {:.3f} + {:.3f}j...'.format(c.real,c.imag))
    sys.stdout.flush()
    for y in range(h):
        for x in range(w):
            z = centre + complex(((x-(w-1)/2.0)*pxw),(((h-1)/2.0-y)*pxw))
            n = 0
            while abs(z) < limit and n < 51:
                try:
                    z = pow(z,a) + c
                    n = n + 1
                except ZeroDivisionError:
                    z = limit
            pixel = int(255 * (0.5 + 0.5*math.cos(math.pi*n/51.0)))
            img.put('#{:02x}{:02x}{:02x}'.format(pixel,pixel,pixel),(x,y))
    sys.stdout.write('OK\n')
    sys.stdout.flush()

def set_z0(event):
    global z0
    z0 = complex((event.x - w/2) * pxw, (h/2 - event.y) * pxw)
    paint()

def set_c(event):
    global c
    c = complex((event.x - w/2) * pxw, (h/2 - event.y) * pxw)
    generate_fractal()
    paint()

def paint():
    canv.delete("all")
    canv.create_rectangle(0, 0, w, h, fill="white")
    canv.create_image((0,0), anchor="nw", image=img, state="normal")
    z = z0
    x = int(w/2 + z.real/pxw)
    y = int(h/2 - z.imag/pxw)
    canv.create_oval(x-3,y-3,x+4,y+4,fill="green")
    for n in range(100):
        p1 = z
        for m in range(M+2):
            if m<=M:
                p2 = pow(z,1.0+((a-1.0)*m)/M)
            else:
                p2 = p2 + c
            x1 = int(w/2 + p1.real / pxw)
            y1 = int(h/2 - p1.imag / pxw)
            x2 = int(w/2 + p2.real / pxw)
            y2 = int(h/2 - p2.imag / pxw)
            canv.create_line(x1, y1, x2, y2, fill="blue")
            p1 = p2
        z = pow(z,a) + c
        x = int(w/2 + z.real/pxw)
        y = int(h/2 - z.imag/pxw)
        canv.create_oval(x-3,y-3,x+4,y+4,fill="red")
        if (abs(z) > 10): break

# Create master Tk widget
master = Tk()
master.title('Julia Set Explorer - by Ted Burke - see http://batchloaf.com')

# Dimensions and resolution
w,h = 800,800
centre = 0 + 0j
pxw = 0.005

# Create Tk canvas widget
canv = Canvas(master, width=w, height=h)
canv.pack()
canv.bind("<B1-Motion>", set_z0)
canv.bind("<Button-1>",  set_z0)
canv.bind("<Button-2>",  set_c)
canv.bind("<Button-3>",  set_c) # this is actually the right click

# Iterating function parameters
a = 2.0
c = 0.325 + 0.285j
z0 = 0 + 0j
limit = 10

# Number of steps in plotting each curve within each orbit
M = 40

# Create image object to store Julia Set image
img = PhotoImage(width=w, height=h)

# Generate initial Julia Set image
generate_fractal()
paint()

# Main Tk event loop
mainloop()

To use the program,

  • Left click anywhere on the image to set the value of z0, which will update the displayed orbit.
  • Clicking and dragging will update z0 and its orbit in real time, which really helps to see what’s going on.
  • Right click anywhere on the image to set c, which will updated the background filled Julia Set image (the black shape). It takes a little time to generate the new image, so be patient!

I’m running Xubuntu 15.04 Linux here, so Python 2.7.9 was already installed by default. I needed to install Tkinter though, which is Python’s standard GUI toolkit. To install Tkinter, I just did

sudo apt-get install python-tk
Posted in Uncategorized | Leave a comment
output
#
# vary_ab.py - Written by Ted Burke, last updated 18-1-2016
#

import subprocess
import math

centre = 0 + 0j   # complex value at centre of image
w,h = 400,400     # image width and height
pxw = 0.02        # pixel width
limit = 4.0       # once z reaches this value, iteration ceases

amin,amax = 2.0,4.0
bmin,bmax = 2.0,4.0

# Create a buffer to store a row of pixel values
row = w*[0]


# Generate images
for frame in range(100):
    # Open file
    filename = 'frame{:03d}'.format(frame)
    pgmfile = file(filename + '.pgm', 'w')
    pgmfile.write('P5\n{} {}\n255\n'.format(w,h))
    for y in range(h):
        for x in range(w):
            a = ((amin+amax)/2.0) + ((amax-amin)/2)*math.cos(2*math.pi * (float(x)/w + frame/100.0))
            b = ((bmin+bmax)/2.0) + ((bmax-bmin)/2)*math.sin(2*math.pi * (float(y)/h + frame/100.0))
            c = centre + complex(((x-(w-1)/2.0)*pxw),(((h-1)/2.0-y)*pxw))
            z = 0
            n = 0
            while abs(z) < limit and n < 51:
                try:
                    z = (pow(z,a) + c)/(pow(z,b) - c)
                    n = n + 1
                except ZeroDivisionError:
                    z = limit
            row[x] = int(255 * (0.5 + 0.5*math.cos(math.pi*n/51.0)))
        pgmfile.write(bytearray(row))

    # Close file
    pgmfile.close()
    subprocess.call('mogrify -format png ' + filename + '.pgm', shell=True)
    subprocess.call('rm ' + filename + '.pgm', shell=True)
    print(filename + '.png written')

#
# renamer.py - Written by Ted Burke, last updated 18-1-2016
#

import subprocess

# Create fractal image files
for n in range(100):
    filename = 'frame{:03d}'.format(n)
    print(filename)    
    subprocess.call('convert ' + filename + '.png -resize 50% ' + 'half_' + filename + '.png', shell=True)

To create an mp4 video using ffmpeg:

ffmpeg -framerate 10 -i half_frame%03d.png -c:v libx264 -pix_fmt yuv420p output.mp4

To create a looping gif using ImageMagick’s convert tool:

convert -loop 0 half_frame*.png output.gif
Posted in Uncategorized | Leave a comment

What I’m working on right now…

I’ve been experimenting with creating fractal images using iterating functions of the form

z_{n+1} = \frac{z_n^a + c}{z_n^b - c}
where
z_0 = 0.

Without getting into all the details right now, within each frame of the animation every pixel represents a different value of the complex variable c. For a given value of c, the number of iterations required for the magnitude of z_n to pass a pre-defined threshold determines the colour of the corresponding pixel. If the sequence exceeds the limit within a small number of iterations, the pixel will be white. Conversely, if the magnitude of z_n remains within the threshold until the maximum number of iterations is reached, the pixel ends up black.

In the following animation, a remains constant (a=5.0) from frame to frame while b is incremented in steps of 0.01 from 0.0 to 4.0. The region of the complex plane shown is -2 < Re{c} < 0 and 0 < Im{c} < 2.

I used the following Python script to generate the individual frames of the animation:

#
# fractal_movie.py - Written by Ted Burke, last updated 7-1-2016
#

import subprocess
import math

centre = 0 + 2j   # complex value at centre of image
w,h = 800,400     # image width and height
pxw = 0.01        # pixel width
limit = 4.0       # once z reaches this value, iteration ceases

# Create a buffer to store a row of pixel values
row = w*[0]

# Create fractal image files
for t in range(401):
    
    a = 5
    b = t/100.0
    
    # Open file
    pgmfilename = 'crpf_a{:1.3f}_b{:1.3f}_pxw{:.3f}_lim{:1.1f}_{:d}x{:d}px.pgm'.format(a, b, pxw, limit, w, h)
    print(pgmfilename)
    pgmfile = file(pgmfilename, 'w')
    pgmfile.write('P5\n{} {}\n255\n'.format(w,h))
    
    # Generate image
    for y in range(h):
        for x in range(w):
            c = centre + complex(((x-(w-1)/2.0)*pxw),(((h-1)/2.0-y)*pxw))
            z = 0
            n = 0
            while abs(z) < limit and n < 51:
                try:
                    z = (pow(z,a) + c)/(pow(z,b) - c)
                    n = n + 1
                except ZeroDivisionError:
                    z = limit
            row[x] = int(255 * (0.5 + 0.5*math.cos(math.pi*n/51.0)))
        pgmfile.write(bytearray(row))
    
    # Close file
    pgmfile.close()
    
    subprocess.call('mogrify -format png ' + pgmfilename, shell=True)
    subprocess.call('rm ' + pgmfilename, shell=True)

I used a second Python script to rename the PNG frames as frame001.png, frame002.png, etc.

#
# renamer.py - Written by Ted Burke, last updated 12-1-2016
#

import subprocess

# Create fractal image files
for t in range(401):
    b = t/100.0
    origfilename = 'crpf_a5.000_b{:1.3f}_pxw0.010_lim4.0_800x400px.png'.format(b)
    newfilename = 'frame{:03d}.png'.format(t)
    print(origfilename + ' ' + newfilename)
    subprocess.call('cp ' + origfilename + ' ' + newfilename, shell=True)

Finally, I used ffmpeg to combine the individual PNG frames into a single video file, as follows:

ffmpeg -framerate 10 -i frame%03d.png -c:v libx264 -pix_fmt yuv420p out.mp4
Posted in Uncategorized | Leave a comment

Simple 2-channel hardware PWM example for the MSP430G2452 microcontroller

//
// 2-channel hardware PWM example for MSP430G2452
// Written by Ted Burke - last updated 2-11-2015
//

#include <msp430.h>
 
int main( void )
{
    // Watchdog timer
    WDTCTL = WDTPW + WDTHOLD; // Disable watchdog timer
    
    // Configure P1.2 (TA0.1) and P1.4 (TA0.2) for PWM output
    // To configure P1.2 as TA0.1: P1DIR.2=1, P1SEL.2=1, P1SEL2.2=0
    // To configure P1.4 as TA0.2: P1DIR.4=1, P1SEL.4=1, P1SEL2.4=1
    P1DIR  = 0b00010100;
    P1SEL  = 0b00010100;
    P1SEL2 = 0b00010000;
    
    // Configure Timer A
    TACTL = TASSEL_2 + ID_0 + MC_1; // Up mode, input divider=1, SMCLK clock
    TACCTL0 = OUTMOD_7;             // TA0.0 in Reset/Set output mode
    TACCTL1 = OUTMOD_7;             // TA0.1 in Reset/Set output mode
    TACCTL2 = OUTMOD_7;             // TA0.2 in Reset/Set output mode
    TACCR0  = 20000;                // Start timer with 20ms period
    
    while(1)
    {
        // P1.2 dim and P1.4 bright for 1 second
        TACCR1 = 2000;
        TACCR2 = 8000;
        __delay_cycles(1000000);
        
        // P1.2 bright and P1.4 dim for 1 second
        TACCR1 = 8000;
        TACCR2 = 2000;
        __delay_cycles(1000000);
    }
    
    return 0;
}
Posted in Uncategorized | 2 Comments

RoboSlam @ Dublin Maker – only two days away!

Looking forward to Dublin Maker on Saturday in Trinity! Entry is free and all are welcome. Our stand will be the one where humans are getting tricked into building a robot army.

RoboSlam

Ted's robot design for Dublin Maker event on Saturday at Trinity College Dublin. Come build one for free! And, it's just 12 Euro (less than the cost to us) if you want to take it home with you. Ted’s robot design for Dublin Maker event on Saturday at Trinity College Dublin. Come build one for free! And, it’s just 12 Euro (less than the cost to us) if you want to take it home with you.

Hello RoboSlammers,

This is just a quick reminder – it is only two days now to Saturday’s big event – Dublin Maker! The weather forecast is looking good so far – “cool and dry with scattered showers” according to met.ie. The event will open at 10am and run until 6pm. It tends to get busy at around midday. A description of some of the main participants at Dublin Maker 2015 is available here. Personally, we can’t wait to see what mechanical wonders Michal Mizsta the “dragon dude”, will have on show this year. And there is an interesting range of exhibits.

Shannon, Damon, Frank, and Ted getting ready today for the RoboSlam Cafe. Shannon, Damon, Frank, and Ted getting ready today…

View original post 239 more words

Posted in Uncategorized | Leave a comment

Generating antiphase PWM signals with the dsPIC30F4011

I frequently receive queries from people who are using a dsPIC microcontroller to control power electronics of some kind, such as in an inverter, a voltage converter, or similar. Many of these queries relate to the generation of different combinations of pulse-width modulation (PWM) signals. In this article, I describe a simple example application in which the dsPIC30F4011 is used to generate two antiphase PWM signals with duty cycle varying between 10% and 45%. By antiphase, I mean that the two signals are 180 degrees (π radians) out of phase with each other. Apart from that, the two signals are identical.

I developed this example in response to a query from gunz which was posted as a comment on another article that I published on this blog a couple of years ago. In that article, I also generated antiphase PWM signals, but I used two of the dsPIC30F4011’s three PWM channels. Specifically, the signals were generated on pins PWM1H and PWM2L. Gunz asked how to generate the signals on the high and low pins of a single PWM channel (e.g. PWM1H and PWM1L), so that’s what I’m demonstrating in this example.

The signal specification is:

  • Two identical PWM signals 180 degrees out of phase.
  • The PWM frequency is 15kHz.
  • The duty cycle varies between 10% and 45% (hence the pulses on the two outputs never overlap).
  • The output pins are PWM1H and PWM1L.

The basic code is as follows:

//
// This dsPIC30F4011 example program generates PWM signals
// on PWM1H and PWM1L. The signals are identical except that
// they are 180 degrees out of phase. Duty cycle can vary
// between 0.1 and 0.45. PWM frequency is constant at 15kHz.
//
// The basic idea is to configure the PWM frequency to double
// what we want (30 kHz), but then only generate a pulse on
// each output every second period. We use the PWM interrupt
// (which is triggered every time the PWM timebase reaches
// the value PTPER and goes back to zero) to switch back and
// forth between PWM1H and PWM1L. i.e. The two pins take it
// in turns to output pulses.
//
// Written by Ted Burke - last updated 11-4-2015
//

#include <xc.h>
#include <libpic30.h>

// Configuration settings
_FOSC(CSW_FSCM_OFF & FRC_PLL16); // Fosc=16x7.5MHz, Fcy=30MHz
_FWDT(WDT_OFF);                  // Watchdog timer off
_FBORPOR(MCLR_DIS);              // Disable reset pin

//
// The PWM interrupt is triggered each time PTMR reaches the same
// value as PTPER and resets to zero. We use this ISR to switch back
// and forth between outputting a pulse on PWM1H and outputting a pulse
// on PWM1L. During each PTMR period, we override one or other pin to
// lock it at 0, suppressing its pulse. Because the OSYNC bit is set,
// PWM override changes do not take effect until the next time PTMR
// resets to zero. Updates therefore don't take effect until the start
// of the the next cycle.
//
void __attribute__((interrupt, auto_psv)) _PWMInterrupt(void)
{
    // Reset PWM interrupt flag
    _PWMIF = 0;
    
    // Alternate between overriding PWM1L and overriding PWM1H.
    // Whichever pin is overridden will be driven low during
    // the next PWM period, suppressing its output pulse.
    if (_POVD1H)
    {
        _POVD1H = 0;
        _POVD1L = 1;
    }
    else
    {
        _POVD1L = 0;
        _POVD1H = 1;
    }
}

int main(void)
{
    // Configure RD0 as a digital output for an indicator LED
    TRISD = 0b1110;
 
    // Configure PWM
    //
    // Note: To provide higher PWM duty cycle resolution, the dsPIC's
    // PDCx unit is only half as long as its PTPER unit. For example,
    // when PDC1 = PTPER, the PWM1 duty cycle is actually only 50%.
    // Furthermore, in this example, since we are alternating back and
    // forth between pulses on two outputs, each output only produces
    // every second pulse. Hence, in this case PTPER is really only
    // half of the full signal period. Yes, it's potentially confusing!
    //
    _PMOD1 = 1;   // Enable PWM channel 1 in independent mode
    _PEN1H = 1;   // Enable PWM1H pin
    _PEN1L = 1;   // Enable PWM1L pin
    _POUT1L = 0;  // When PWM1L is overriden, set it low
    _POUT1H = 0;  // When PWM1H is overriden, set it low
    _POVD1L = 1;  // Initially, enable override on PWM1L
    _POVD1H = 0;  // Initially, disable override on PWM1H
    _OSYNC = 1;   // Synchronise PWM override changes with PTMR reset
    _PWMIE = 1;   // Enable PWM interrupt
    _PTCKPS = 0;  // prescale=1:64 (0=1:1, 1=1:4, 2=1:16, 3=1:64)
    PTPER = 999;  // Set PWM frequency to 30 kHz
    PDC1 = 400;   // Set duty cycle to 10%
    _PTEN = 1;    // Enable PWM time base
    
    // Now just blink an LED while the PWM ISR does the heavy lifting
    while(1)
    {
        _LATD0 = 1;          // Turn on LED on RD0
        __delay32(15000000); // 0.5 second delay
        _LATD0 = 0;          // Turn off LED on RD0
        __delay32(15000000); // 0.5 second delay
    }

    return 0;
}

I compiled the above code using Microchip’s free XC16 C compiler. This is my build script, which produces an binary file called “a.hex”:

xc16-gcc main.c -mcpu=30F4011 -Wl,--script=p30F4011.gld
if errorlevel 0 xc16-bin2hex a.out

I downloaded the file “a.hex” onto the dsPIC30F4011 using a PICkit 2 USB programmer. I don’t use MPLAB at all, so I just used the PICkit 2 software application (V 2.61 can be downloaded here).

The breadboard circuit I used for testing is shown below. In addition to the voltage supply connections I included in my circuit, I highly recommend connecting pin 20 (VSS) and pin 21 (VDD) to the voltage supply rails. In the past, I have tended not to bother connecting all voltage supply pins, but I have recently run into some very strange problems that it took me a long time to realise could be completely solved by connecting more voltage supply pins. Specifically, I had a lot of problems with mysterious resetting of the dsPIC when using PWM interrupts in parallel with mainline code (such as the LED flashing code I put in the main function in this example). I therefore highly recommend connecting all of the dsPIC’s voltage supply pins, even if it seems redundant.

An indicator LED (with 220\Omega; current-limiting resistor in series) is connected to RD0 (pin 23). The three wires which extend beyond the lower edge of the image are the oscilloscope connections – green is ground and the two blue wires connect PWM1H and PWM1L to different channels of the scope.

20150410_181128

This is how the circuit appeared once the code shown above was running. The LED on RD0 simply blinks on and off once a second.

I displayed the signals from PWM1H and PWM1L on the two channels of an oscilloscope. As shown below, the two signals are identical but 180 degrees out of phase.

Antiphase PWM signals displayed on oscilloscope.

The following modified version of the previous example shows how the PWM ISR can be used to update the duty cycle as well as alternating the generated pulses between the two outputs. Apart from comments, the only change from the previous example is the addition of lines 62-68 to the PWM ISR (the _PWMInterrupt function) which modulate the PWM duty cycle.

//
// This dsPIC30F4011 example program generates PWM signals
// on PWM1H and PWM1L. The signals are identical except that
// they are 180 degrees out of phase. This is a slightly
// modified version of the example in which duty cycle is
// modulated in a sawtooth pattern, increasing in small steps
// from 0.1 all the wat to 0.45, then resetting to 0.1 over
// and over again. PWM frequency is constant at 15kHz.
//
// The basic idea is to configure the PWM frequency to double
// what we want (30 kHz), but then only generate a pulse on
// each output every second period. We use the PWM interrupt
// (which is triggered every time the PWM timebase reaches
// the value PTPER and goes back to zero) to switch back and
// forth between PWM1H and PWM1L. i.e. The two pins take it
// in turns to output pulses.
//
// Written by Ted Burke - last updated 11-4-2015
//

#include <xc.h>
#include <libpic30.h>

// Configuration settings
_FOSC(CSW_FSCM_OFF & FRC_PLL16); // Fosc=16x7.5MHz, Fcy=30MHz
_FWDT(WDT_OFF);                  // Watchdog timer off
_FBORPOR(MCLR_DIS);              // Disable reset pin

//
// The PWM interrupt is triggered each time PTMR reaches the same
// value as PTPER and resets to zero. We use this ISR to switch back
// and forth between outputting a pulse on PWM1H and outputting a pulse
// on PWM1L. During each PTMR period, we override one or other pin to
// lock it at 0, suppressing its pulse. Because the OSYNC bit is set,
// PWM override changes do not take effect until the next time PTMR
// resets to zero. Updates therefore don't take effect until the start
// of the the next cycle.
//
void __attribute__((interrupt, auto_psv)) _PWMInterrupt(void)
{
    // Reset PWM interrupt flag
    _PWMIF = 0;
    
    // Alternate between overriding PWM1L and overriding PWM1H.
    // Whichever pin is overridden will be driven low during
    // the next PWM period, suppressing its output pulse.
    if (_POVD1H)
    {
        _POVD1H = 0;
        _POVD1L = 1;
    }
    else
    {
        _POVD1L = 0;
        _POVD1H = 1;
    }

    // This addition to the PWM ISR modulates the duty cycle in
    // a sawtooth pattern. Every 100th time the ISR runs, the
    // duty cycle is increased by incrementing PDC1. When the
    // duty cycle reaches 45%, it is reset to 10%.
    static int n=0;
    n = (n + 1) % 100;
    if (n==0)
    {
        PDC1 = PDC1 + 1;
        if (PDC1 > (0.45 * 4 * PTPER)) PDC1 = 0.1 * 4 * PTPER;
    }
}

int main(void)
{
    // Configure RD0 as a digital output for an indicator LED
    TRISD = 0b1110;
 
    // Configure PWM
    //
    // Note: To provide higher PWM duty cycle resolution, the dsPIC's
    // PDCx unit is only half as long as its PTPER unit. For example,
    // when PDC1 = PTPER, the PWM1 duty cycle is actually only 50%.
    // Furthermore, in this example, since we are alternating back and
    // forth between pulses on two outputs, each output only produces
    // every second pulse. Hence, in this case PTPER is really only
    // half of the full signal period. Yes, it's potentially confusing!
    //
    _PMOD1 = 1;   // Enable PWM channel 1 in independent mode
    _PEN1H = 1;   // Enable PWM1H pin
    _PEN1L = 1;   // Enable PWM1L pin
    _POUT1L = 0;  // When PWM1L is overriden, set it low
    _POUT1H = 0;  // When PWM1H is overriden, set it low
    _POVD1L = 1;  // Initially, enable override on PWM1L
    _POVD1H = 0;  // Initially, disable override on PWM1H
    _OSYNC = 1;   // Synchronise PWM override changes with PTMR reset
    _PWMIE = 1;   // Enable PWM interrupt
    _PTCKPS = 0;  // prescale=1:64 (0=1:1, 1=1:4, 2=1:16, 3=1:64)
    PTPER = 999;  // Set PWM frequency to 30 kHz
    PDC1 = 400;   // Set duty cycle to 10%
    _PTEN = 1;    // Enable PWM time base
    
    // Now just blink an LED while the PWM ISR does the heavy lifting
    while(1)
    {
        _LATD0 = 1;          // Turn on LED on RD0
        __delay32(15000000); // 0.5 second delay
        _LATD0 = 0;          // Turn off LED on RD0
        __delay32(15000000); // 0.5 second delay
    }

    return 0;
}

Here’s a short video showing the sawtooth modulation of the PWM duty cycle displayed on the oscilloscope.

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