Heaviside Step Function

I would like to plot this function but I’m not sure how to implement it.
It is a Heaviside step function, and I assume that the theta sign means that the sinusoid applies only to 27π < t < 31π , and otherwise the function would be zero.
Any hint?

1 Like

Where did you get the screenshot from?

This is the whole parametric function.

I must be wrong with the assumption of the probably meaning of the UnitStep function. because the butterfly plot (link above), plots t from 0 to 27*PI.
So I can’t write:

float dt =  PI/3000; // stepping
  for (float t = 0; t <= 27*PI; t += dt) {
    if (t < 31*PI && t > 27*PI) x = sin(9*t+14/3)+247; 
    else x = 0; 

because it only increments until 27*PI

That’s wild. Are you trying to implement the entire function in code?

Not that wild, It´s just complex. But you can start with simple ones
Let’s take a random simple parametric curve function of their site like this one.
The code in processing is:

float dt; // delta stepper used in for loop. has to be small enough to give a good resolution
float x, px, y, py;  // actual and previous x/y values
int sf; // scale factor

void setup() {
  size(800, 800);
  sf = 300; 
  dt = PI/1500;
  translate(width/2, height/2);
  for (float t = dt; t < TWO_PI; t += dt) {
    x = sin(10*t);
    y = sin(8*t); 
    line(sf*px, -sf*py, sf*x, -sf*y);
    px = x;
    py = y;

The complex ones like the Einstein figure work the same way.
But they use the step function to omit/skip/hide certain parts of the curves. End that is the part with what I’m struggling with.
The WolframAlpha site is really interesting. You can abstract functions from images etc. But it’s paid. An introduction you can find here and here.
There you can find fi the image below made with solely functions.


1 Like

@noel, I’m not sure about the theta issue specifically, but have you looked at @quark’s Jasmine library?


It might be possible for you to feed in your equation components in something more like their natural equation form, rather than hand-translating them into Java functions.

That said, I haven’t tried a project like this with Jasmine before – it just occurred to me that it might be possible and labor saving, especially with so many expressions in things like person curves.


Interesting problem :smile:

Some points

  • the equation you present is invalid because there is an imbalance in the brackets, I assume that there is an opening bracket before the sin
  • the function as presented will return 0 outside the range 27\pi to 31\pi
  • inside the range 27\pi to 31\pi it will return a sinusoidal value between 246 and 248

The solution is to create another function for \theta(v) which in the code below is called h(v)

void setup() {
  for (float t = 26.9*PI; t < 31.1*PI; t += 0.1) {
    println(t/PI, func(t));

float func(float t) {
  float result = sin(9*t + 14f/3f) + 247;
  result *= h(31 * PI - t);
  result *= h(t - 27 * PI);
  return result;

// Calculate the Heavyside step function H((t + n*Pi) 
// Return 0 if t + n*Pi <0
// Return 1 if t + n*Pi >0
// otherwise return 0.5
float h(float v) {
  return v < 0 ? 0 : v > 0 ? 1 : 0.5;

Hi @quark .
Thank you for looking into this and providing this function.
Unfortunately I have still difficulties in implementing it the way as done on the wolframalpha site.
It´s not that I want to plot the existing functions on their site, but rather trying to discover how they derive those functions from an image. I do not want to pay monthly to use their programs and learn yet another language (the wolfram language); but by searching throughout their blogs I hope to discover something.
I would really appreciate it, if you could have a look on one of their blogs. Scrolling down one third of this blog-page, you will find a plot function of the twitter bird. As you see below I coded this shape in a para-metrical way. How can I use your step function in this code?
Any hint about how they derive those functions will be a great help.
Thanks once again.

float[][] cc = {{210,168},{400,148},{210,280},{70,185},{130,215},{180,216},{130,110},{185,135},{292,-75},{430,30},{493,-35},{410,10},{450,-75}}; // circle center
float[][] ca = {{-.15*PI,.5*PI},{.66*PI,1.56*PI},{-.38*PI,0*PI},{0*PI,.2*PI},{-.07*PI,.06*PI},{-.54*PI,-.08*PI}, {-.18*PI,-.08*PI},
                {-.67*PI,-.23*PI},{-.29*PI,-.02*PI},{.18*PI,.3*PI},{0.02*PI,.12*PI},{.22*PI,.35*PI},{.08*PI,.16*PI}}; // circle arc {start,stop}
int[] radius = {300,105,110,250,110,110,110,110,255,148,160,150,177};

void setup() {
  size(620, 520);
  stroke(50, 50, 255);
  float step = 2*PI/60;
  float px = 0, py = 0;
  for (int i = 0; i < 13; i++) {
    translate(cc[i][0], cc[i][1]);
    int j = 0;
    for (float t = ca[i][0]; t < ca[i][1]+step; t += step) {
      float x = radius[i]*sin(t);
      float y = radius[i]*cos(t);
      if (j == 0) { px = x; py = y; } // avoiding -center to start line- in first for loop
      line(x, y, px, py);
      px = x; py = y; j++;
  // fillShape();
Not relevant fillShape function
import java.awt.Point;
import java.util.Queue;
import java.util.LinkedList;
PImage img;

void fillShape() {
  img = get();
  Queue<Point> queue = new LinkedList<Point>();
  queue.add(new Point(width/2, height/2));
  while (!queue.isEmpty()) {
    Point p = queue.remove();
    if (check(p.x, p.y)) {     
      queue.add(new Point(p.x, p.y-1)); 
      queue.add(new Point(p.x, p.y+1)); 
      queue.add(new Point(p.x-1, p.y)); 
      queue.add(new Point(p.x+1, p.y));
  image(img, 0, 0, width, height);

boolean check(int x, int y) {
  if (x < 0 || y < 0 || y >= img.height || x >= img.width) return false;
  int pp = img.pixels[x+(y*img.width)];
  if (pp != color(255)) return false; 
  img.pixels[x + (y * img.width)] = color(0, 167, 230); 
  return true;

1 Like

They seem to use a wide range to techniques depending on the image type and the maths is hidden inside the Wolfram language.

The HSF is used to limit the parts of individual curves to include in the calculation so is useful if you are trying to describe the shape using a single pair of parametric curves i.e. x=f(t) and y=g(t) which is what they have in the Einstein face drawing.

In your code you are creating the twitter outline using a series of arcs and using the ca array to control the portion of the circle to draw so you don’t have a need for the HSF.

The approach you take from here depends on whether you want to

  1. create a single pair of parametric equations . x=f(t) and y=g(t) to describe the whole contour
  2. describe the contour as a collection of circular arcs

If you want (1) then you need to modify the data so that the contour is roughly centred about the coordinates [0,0]. If the contour is offset from the origin then most values of the parametric parameter t would cause the HSV to output 0 which is hugely wasteful.

If you want (2) I suggest that you create a class that represents an arc it would have fields for the centre, radius and angle limits. Then have an array of arcs, this is better than having separate arrays for centre, angle limits and radii.

Personally I would go for (2) because it is easy to create a small sketch that takes an image and with some user input calculate the arc values need for the contour.

When drawing the arcs instead of translating the graphics matrix you calculate the arc cartesian values with
x = cx + radius * sin(t) and
y = cy + radius * cos(t)
it should be possible to generate the parametric equations for method (1), though I would need to think more on the actual algorithm.


Thanks. Great hint. From the center! That immediately reminded me of Ptolemy (proving that that the earth is the center of the universe). While I took the second approach, WolframAlpha likely uses the first, using Fourier Transform in their math. So I did a search and found this!

See the Pen Sketchable Fourier Transform by marl0ny (@marl0ny) on CodePen.

Or using ready shapes:

See the Pen Drawing Paths with Harmonics by Mitch Anderson (@tmanderson) on CodePen.

Edit: I just found the series Daniel Shiffman made about this subject.
Plenty of material to study nowl
Starting here.

1 Like

I really enjoyed studying this. Daniel really is a great teacher, Everything is in P5.js and I couldn’t find a ready code for the ‘fourier complex’ sketch example in P5.java, so I will leave it here. It also includes a way of evenly distributing points on the curves. By using parametric functions, the huge data file is no longer necessary.
So did I come closer on discovering how Wolfram derive those functions from an image? Nope.
But paying more attention, I noticed that their functions use as many quotients as I use floats, in my code, and that is what I’m looking after; plotting images with the least amount of data.
So my next project is drawing free hand curves, and transform them in bezier curves.
If anyone can give some guide lines on how to start doing this I would appreciate it.

//  Code based on teachings from Daniel Sauter https://www.youtube.com/watch?v=7_vKzcgpfvU&t=1314s

float[][] cc = {{210, 168}, {430, 30}, {493, -35}, {410, 10}, {450, -75}, {400, 148}, {292, -75}, {185, 135}, {130, 110}, {180, 216}, {130, 215}, {210, 280}, {70, 185}}; // circle center
float[][] ca = {{-.15, .5}, {.18, .3}, {0.02, .12}, {.22, .35}, {.08, .16}, {.66, 1.56}, {-.29, -.02}, {-.67, -.23}, {-.18, -.08}, {-.54, -.08}, {-.07, .06}, {-.38, 0}, {0, .2}}; // circle arc {start,stop}
int[] radius = {300, 148, 160, 150, 177, 105, 255, 110, 110, 110, 110, 110, 250};
float[] bird_x = new float[0];
float[] bird_y = new float[0];
float[] swap_x = new float[0];
float[] swap_y = new float[0];
int count = 0, tel;
float diameter, amplitude, time, interval, frequency, phase;
PVector pos, prev;
ArrayList<PVector> path, drawing;
ArrayList<Complex> x, fourier;
boolean allow = true;

void setup() {
  size(800, 600);
  pos = new PVector();
  prev = new PVector();
  x = new ArrayList<Complex>();
  fourier = new ArrayList<Complex>();
  path = new ArrayList<PVector>();
  float step = .037;
  float px = 0, py = 0;
  int tel = 0;
  for (int i = 0; i < 13; i++) {
    translate(cc[i][0], cc[i][1]);
    int j = 0;
    swap_x = new float[0];
    swap_y = new float[0];
    for (float t = ca[i][0]*PI; t < ca[i][1]*PI+1.7*step; t += step/radius[i]*200) {
      float x = radius[i]*sin(t);
      float y = radius[i]*cos(t);
      swap_x = append(swap_x, x+cc[i][0]);
      swap_y = append(swap_y, y+cc[i][1]);
    if (i != 0 && i % 2 == 0) {
      swap_x = reverse(swap_x); 
      swap_y = reverse(swap_y);
    bird_x = concat(bird_x, swap_x);
    bird_y = concat(bird_y, swap_y);
  PVector[] vectors = new PVector[bird_x.length]; 
  for (int i = 0; i < bird_x.length; i++) {
    vectors[i] = new PVector();
    vectors[i].x = bird_x[i];
    vectors[i].y = bird_y[i];
  drawing = new ArrayList<PVector>();
  for (int i = 0; i < vectors.length; i++) {
  for (int i = 0; i < drawing.size()-1; i++) {
    PVector point = drawing.get(i);
    Complex c = new Complex(point.x, point.y);
  fourier = dft(x);

void draw() {
  PVector vertex = drawEpicycles(width/2, height/2, fourier);
  stroke(0, 0, 255);
  for (int i = path.size()-1; i > 0; i--) {
    PVector p = path.get(i);
    vertex(p.x, p.y);
  time += interval;
  if (time > TWO_PI) {
    time = 0;

PVector drawEpicycles(float x, float y, ArrayList<Complex> fourier) {
  pos.x = x;
  pos.y = y;
  interval = TWO_PI/fourier.size();
  for (int i = 1; i < fourier.size(); i++) {
    prev.x = pos.x;
    prev.y = pos.y;
    Complex epicycle = fourier.get(i);
    frequency = epicycle.freq;
    amplitude = epicycle.amp;
    phase = epicycle.phase;
    diameter = amplitude*3;
    float theta = frequency * time + phase;
    pos.x += amplitude * cos(theta);
    pos.y += amplitude * sin(theta); 
    stroke(0, 150, 0);
    if (i == 1) diameter = 400;
    ellipse(prev.x, prev.y, diameter, diameter);
    line(prev.x, prev.y, pos.x, pos.y);
  return new PVector(pos.x, pos.y);

ArrayList<Complex> dft(ArrayList<Complex> x) {
  int N = x.size();
  ArrayList<Complex> X = new ArrayList<Complex>(N);
  for (int k = 0; k < N; k++) {
    Complex sum = new Complex(0, 0);
    for (int n = 0; n < N; n++) {
      float phi = (TWO_PI * k * n) / N;
      Complex c = new Complex(cos(phi), -sin(phi));
      sum = sum.add(x.get(n).mult(c));
    sum = sum.div(N);
    float freq = k;
    float amp = sum.mag();
    float phase = sum.heading();
    X.add(new Complex(sum.re, sum.im, freq, amp, phase));
  return X;

void sort(ArrayList<Complex> c) {
  int n = c.size();
  for (int i = 0; i < n-1; i++) {
    int temp = i;
    for (int j = i+1; j < n; j++) {
      if (c.get(j).amp > c.get(temp).amp) temp = j;
    Complex temp2 = c.get(temp);
    c.set(temp, c.get(i));
    c.set(i, temp2);

class Complex {
  float re;
  float im;
  float freq;
  float amp;
  float phase;

  Complex(float r, float i) {
    re = r;
    im = i;

  Complex(float r, float i, float f, float a, float p) {
    re = r;
    im = i;
    freq = f;
    amp = a;
    phase = p;

  Complex mult(Complex c_obj) {
    float rea = re * c_obj.re - im * c_obj.im;
    float ima = re * c_obj.im + im * c_obj.re;
    return new Complex(rea, ima);

  Complex add(Complex c_obj) {
    return new Complex(re + c_obj.re, im + c_obj.im);

  Complex div(float divisor) {
    float r = re / divisor;
    float i = im / divisor;
    return new Complex(r, i);

  float heading() {
    return atan2(im, re);

  float mag() {
    return sqrt(re * re + im * im);


In addition to Dan’s tutorial, I also enjoy this Alex Miller illustrated explainer:

It isn’t in Processing / p5.js – but the discussion with illustrations is still helpful.

1 Like

The Wolfram people drawing parametric equations used complex numbers!! Now that would be fun.

What do you mean by that?

I meant the challenge to learn why they use complex numbers and recreate the techniques would be fun.

My Steganos library came about after watching the film Along Came a Spider, it was the first time I had heard of steganography (hiding in plain sight) and after some research I thought it would be fun to implement it. So Along Came Steganos and although I have never used the library for real myself it doesn’t spoil the satisfaction of creating it.

1 Like