The FFT is a simple function, that takes as input a periodic signal from one domain and outputs a delta spike in another.


In [12]:
%pylab inline
import pims
from __future__ import print_function
from __future__ import division
N = 512
phs = linspace(0.2* pi, 7.8 * 2 * pi, N)
x = sin(phs)
plot(x)
title('A Beautiful Periodic Function');


Populating the interactive namespace from numpy and matplotlib

We start with a beauitful periodic function... perform a bunch of multiplies and adds... they look and mean something like this: $$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi k n / N}$$ Which is the Discrete Fourier Transform written using math.

$$X_k = \sum_{n=0}^{N-1} x_n \cdot [\cos(-2 \pi k n / N) + i \sin(-2 \pi k n / N)]$$$$X_k = \sum_{n=0}^{N-1} x_n \cdot [\cos(2 \pi k n / N) - i \sin(2 \pi k n / N)]$$

In [13]:
X = fft.rfft(x)
plot(abs(X))
title('The real part of our FFT on the Beautiful Periodic Function.');


Remember that, even though this seems like a tiny function or a rather small mathematical operation, it still looks like this written in C code.

Source: Paul Bourke</A> / This computes an in-place complex-to-complex FFT x and y are the real and imaginary arrays of 2^m points. dir = 1 gives forward transform dir = -1 gives reverse transform / short FFT(short int dir,long m,double x,double y) { long n,i,i1,j,k,i2,l,l1,l2; double c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++) 
      n *= 2;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }

   return(TRUE);
}

Remember one more thing, the C code still is quite likely to travel through a single core CPU, designed using the x86 Instruction Set Architecture, that for an old model of processor design might look something like this. (Source: 8086 Architecture )


In [21]:
j = pims.ImageSequence('imgs/*.jpg')
%matplotlib inline
j[0]


Out[21]:

We know exactly how to talk to the batch of protons and electrons which are represented by the above Register Transfer Level Diagram. We use these simple instructions:

AAA, AAD, AAM, AAS, ADC, ADD, AND, CALL, CNW, CLC, CLD, CLI, CMC, CMP, CMPSB, CMPSW, CWD, DAA, DAS, DEC, DIV, ESC, HLT, IDIV, IMUL,IN,INC,INT,INTO, IRET, Jcc, JCXZ, JMP, LAHF, LDS, LEA, LES, LOCK, LODSB, LODSW, LOOP/LOOPX, MOV, MOVSB, MOVSW, MUL, NEG, NOP, NOT, OR, OUT, POP, POPF, PUSH, PUSHF, RCL, RCR, REPXX, RET, RETN, RETF, ROL, ROR, SAHF, SAL, SAR, SBB, SCASB, SCASW, SHL, SHR, STC, STD, STI, STOSB, STOSW, SUB, TEST, WAIT, XCHG, XLAT, XOR.

But do not worry about our Instruction Set Architecture. Python handles everything. Just remember, every time something simple happens, like one tiny multiplication of two sixteen bit numbers, this much work is happening.


In [25]:
%%html
<div> Input is Blue output is Yellow. All other colors are logic gates (AND, OR, NOT)</div>
<div id="d3-example"></div>
<style>
.node {stroke: #fff; stroke-width: 1.5px;}
.link {stroke: #999; stroke-opacity: .6;}
</style>


Input is Blue output is Yellow. All other colors are logic gates (AND, OR, NOT)

In [27]:
%%javascript
// We load the d3.js library from the Web.
require.config({paths: {d3: "http://d3js.org/d3.v3.min"}});
require(["d3"], function(d3) {
    // The code in this block is executed when the 
    // d3.js library has been loaded.
    
    // First, we specify the size of the canvas containing
    // the visualization (size of the <div> element).
    var width = 600,
        height = 600;

    // We create a color scale.
    var color = d3.scale.category10();

    // We create a force-directed dynamic graph layout.
    var force = d3.layout.force()
        .charge(-20)
        .linkDistance(30)
        .size([width, height]);

    // In the <div> element, we create a <svg> graphic
    // that will contain our interactive visualization.
    var svg = d3.select("#d3-example").select("svg")
    if (svg.empty()) {
        svg = d3.select("#d3-example").append("svg")
                    .attr("width", width)
                    .attr("height", height);
    }
        
    // We load the JSON file.
    d3.json("sample_circuit/sixteenbitmultiplier.json", function(error, graph) {
        // In this block, the file has been loaded
        // and the 'graph' object contains our graph.
        
        // We load the nodes and links in the force-directed
        // graph.
        force.nodes(graph.nodes)
            .links(graph.links)
            .start();

        // We create a <line> SVG element for each link
        // in the graph.
        var link = svg.selectAll(".link")
            .data(graph.links)
            .enter().append("line")
            .attr("class", "link");

        // We create a <circle> SVG element for each node
        // in the graph, and we specify a few attributes.
        var node = svg.selectAll(".node")
            .data(graph.nodes)
            .enter().append("circle")
            .attr("class", "node")
            .attr("r", 4)
            //.attr("r", function(d){return (d.id)*2} )  // radius
            .style("fill", function(d) {
                // The node color depends on the club.
                return color(d.gate); 
            })
            .call(force.drag);

        // We bind the positions of the SVG elements
        // to the positions of the dynamic force-directed graph,
        // at each time step.
        force.on("tick", function() {
            link.attr("x1", function(d) { return d.source.x; })
                .attr("y1", function(d) { return d.source.y; })
                .attr("x2", function(d) { return d.target.x; })
                .attr("y2", function(d) { return d.target.y; });

            node.attr("cx", function(d) { return d.x; })
                .attr("cy", function(d) { return d.y; });
        });
    });
});


The FFT is acting on our single signal. What happens when we use it on a different type of signal? Like an image? Remember, these bits arrive in parallel in the time domain. This is different from our original signal, which handled a beautiful periodic function which arrived sequentially in the time domain. Go to next notebook