Random noise textures on the GPU in WebGL

Procedural rendering routines often need pseudorandom numbers. For graphics rendering, we usually prefer correct-ish and fast code over properly "correct" code. My approach to random number generation is usually to intermittently seed a weak pseudo-random number generator, like a linear-feedback shift register, from a stronger source like Javascript's Math.random(). 

On the GPU, things are not so simple: we want a pseudorandom algorithm that is fast, local, and parallel. For compatibility, we should also have this RNG use an 8-bit RGBA texture for its internal state. Bitops are a bit tricky in WebGL shaders (though not impossible), so we'll use a linear congruential generator.

Visualizing z-transforms in the the complex plane

These notes (pdf) , along with the script plotting_helper.py , are available as an ipython notebook on Github .

Polar and rectangular coordinates

Recall the rectangular $z=x+ iy$ and polar $z=r e^{i\theta}$ forms of a complex number.

Coloring the points in the complex plane

The two maps below are colored representations of the magnitude (left) and phase (right) in polar coordinates. Observe that the phase rotates $2\pi$ around zero. To connect this to the analysis of poles and zeros in our z-domain transfer function, we can think of these maps as coloring the identity map. For example, it could be a z-domain transfer function of the form: $$ U(z) = z $$ This map has has one zero (at zero), and one pole (at infinity).

The way that phase behaves near 0 and $\infty$ in the identity map resembles the behavior of phase around zeros and poles in more complicated s-domain and z-domain transfer functions. In the identity map, going around the origin counterclockwise adds $2\pi$ to the phase. This is equivalent to encircling infinity clockwise . In a transfer function, each zero (pole) can be viewed as a shifted (for poles: and inverted) version of the identity map.

This allows us to determine the # of poles - # zeros by taking a closed contour integral in the complex plane, which can help us assess the stability of a sytem. For example, in continuous time we need there to be no poles in the right half plane. In discrete time, we need to check that there are no poles outside the unit circle.

This follows from Cauchy's argument principle and is important intution for understanding the Nyquist stability criterion in both the continuous (s-domain) and discrete (z-domain). However, we do not need to go through the proof of the argument principle to understand it intuitively.

Maps in the complex plane

By coloring points in the complex plane, we can visualize functions from $\mathbb C\to\mathbb C$

A shift by 1 (zero at -1)

Note that the $2\pi$ phase rotation around the zero, located at -1

Squaring (two zeros at zero)

Note that the phase rotation doubles, and is now $4\pi$ rather than $2\pi$ around zero.

Inverting (pole at 0)

Note that the phase rotation has reversed.

Let's make an unstable conjugate pair z-domain transfer function

$$ U(z) = \frac 1 {(1 + \alpha e^{i\pi/4} z^{-1})(1 + \alpha e^{-i\pi/4} z^{-1})},\,\, \alpha < 1 $$

Phase rotation vanishes around regions that contain equal # zeros and poles

Poles and zeros cancelling inside unit circle: Bode and Nyquist plots

Phase and magnitude together

We can plot the output of a function as a shaded color, where the hue represents the phase and the brightness the absolute magnitude. (We'll use this in a bit to visualize transormation in the complex plan in a single plot)

Just for fun: Fractals as iterated polynomial maps

Just for fun: the iterated map $z_n = z_n^2 + c$ can generate Julia set fractals. Here's an example for $c = -0.81 + 0.05i$, iterated 7 times.

Approximate conversions from Laplace to Z domain

Forward

$$ z = e^{sT} \approx 1 + sT \Rightarrow sT \approx z-1 $$

Backward

$$ 1/z = e^{-sT} \approx 1 - sT \Rightarrow sT \approx \frac{z-1}{z} $$

Tustin/Bilinear

$$ z = \frac {e^{\tfrac 1 2 sT}}{e^{-\tfrac 1 2 sT}} \approx \frac {1 + \tfrac 1 2 sT} {1 - \tfrac 1 2 sT} \Rightarrow sT \approx 2\frac {z-1}{z + 1} $$

Applying each of these transformations in the complex plane...

Applied to a Laplace domain transfer function w. two conjugate pair poles

• All transfomations introduce some error
• Forward Euler can be unstable
• Backward Euler
  • Maps the point at $\infty$ to 0
  • Preserves stability
  • Can distort frequences and phase response
• Tustin transformation
  • Maps the point at $\infty$ to -1
  • Preserves stability
  • Can distort frequences and phase responses
• Transforming poles and zeros
  • Poorly approximates the linear time-invariant frequency response

Comparing frequency domain responses

Note that, even through the system attained via forward Euler is unstable, it still appears to have a reasonable frequency response. This is because we cannot tell whether the poles are inside or outside the unit circle from the frequency response alone, and highlights the importance of checking the stability before interpreting freqeuency reponse plots.

( In these plots, we use an arbitrary (normalized) frequency, and set the sampling rate of the discrete system to be twice this unit frequency (T=0.5) )

Extra example 1

$$ \frac{10^{-4} z^3}{(z-1)(1.15z^2-2.05z+1)} $$

Extra example 2

$$ \frac 1 {z^2(z-1)} $$