Andrew Que Sites list Photos
Projects Contact
Main

February 04, 2017

AQ Graph Release

   After a fair bit of cleanup work, my Javascript graphing library AQ Graph is ready for the world.
Graph demo   The design of AQ Graph began more than 6 months ago after an experiment at work with HTML 5 canvases.  I had started the basics for a simple line graph but ended up not needing any graphing functionality on the project I was using.  Then some months ago I started working on some FFT code where I needed a Javascript graphing library.  In the past I have used my PHP library X/Y Plot, but I wanted someone that ran client-side rather than server-side.  So I decided to look into what I was going to use for work but never pursued.  The results are a very fast X/Y graphing library that work with large data sets.  As part of a (as yet unreleased) project called WaveLab, I finished the base implementation of the graphing library.  This project required two kinds of graphs: an X/Y plot, and a bar graph.  They were quite messy so I decided that before finishing WaveLab, I would finish the graphing library.  The finishing process involved taking the project specialized library and making it general purpose. 
   I kept the high speed of the X/Y graph and worked on cleaning up the division and labeling function, making them common between the X/Y and bar graph types.  The bar graphs needed the most work.  While fairly basic in terms of options, they are fully functional.  I also made them capable of running as stand-alone and AMD libraries.

January 05, 2017

FFT for LibreOffice Calc

A few days ago I posted a ooBasic script that can be used in LibreOffice to do a Discrete Fourier Transform (DFT) on a set of cells. While it works, it is slow. This is because the a raw implementation of a discrete Fourier transform runs at a complexity of O( n2 ). That is, the number of operations required to compute is the square of the number of data points. Since this grows exponentially, slightly larger data sets take significantly longer to compute.

Naturally, the solution is to switch to a Fast Fourier Transform (FFT). These can generally reduce the complexity down to O( n log2 n ). I knew this when implementing a DFT, but implementing an FFT is a bit more complected. ooBasic kinda...well, pretty much sucks as a programing language. So to do an FFT implementation I wanted an algorithm fairly straight forward to translate. I surveyed several libraries but ultimately decided on Free small FFT. It included a fairly clean radix-2 Cooley–Tukey FFT algorithm but also had Bluestein's algorithm for data not in powers-of-two. The complexity looks to be something around:

I am judging this based on the fact it does three radix-2 transforms on a four times the rounded up power-of-two.

Why does anyone need an FFT in a spreadsheet? I’m not the only person who has searched for such a function, but the primary reason I want one is to do a spectrum plot of a sampled signal. It doesn’t happen often that I need one, but when I want it, I want it in the spreadsheet and would rather not have to go to an external tool.

Today I finished making a simple project page to host the script. Installing it in LibreOffice isn’t easy, but the steps outline in the manual section seem to work. They don’t make it easy, but should anyone else ever want LibreOffice Calc to do an FFT, this should work.

December 27, 2016

LibreOffice Calc Discrete Fourier Transform

      I've been tired of not having one of these, so today I implemented a macro for doing a Discrete Fourier Transform (DFT) in LibreOffice Calc.  I implemented a DFT rather than an Fast Fourier Transform (FFT) because I wanted to be able to use an arbitrary data size, and because the algorithm is pretty easy to implement.  I cannot remember the last time I did any coding in Basic, but I still don't like it.  I find arrays irritating when implemented with the offset of 1 rather than 0.  Guess I'm too use to C.
' Discrete Fourier Transform (DFT).
' The complexity of this function is O( n^2 ) as it is not a Fast Fourier
' Transform (FFT).  However, unlike as FFT this function works an any input size.
Function dft_guts( inputData As Variant, direction As Double, isScale As boolean )
 Dim OFFSET As Integer
 OFFSET = LBound( inputData, 1 )
   
 Dim NUMBER As Integer
 NUMBER = UBound( inputData, 1 ) - OFFSET
 
 ' If we are only given a single column of data, assume that holds the real
 ' values, and the imaginary values are all 0.  Make a new array with the
 ' additional imaginary values in it.
 If UBound( inputData, 2 ) < 2 Then
   Dim newArray( LBound( inputData, 1 ) To UBound( inputData, 1 ), LBound( inputData, 1 ) To UBound( inputData, 1 ) ) As Double
   
   ' Copy data into new array.
   Dim index As Integer
   For index = LBound( inputData, 1 ) To UBound( inputData, 1 )
     newArray( index, 1 ) = inputData( index, 1 )
     newArray( index, 2 ) = 0.0
   Next
   
   ' Switch to new array.
   inputData = newArray
 EndIf

 ' Storage for results.
 Dim outputData( 0 To NUMBER, 1 ) As Double
 
 Dim multiplier As Double
 multiplier = direction * 2 * Pi
 
 ' If scaled, the scale is the number of samples.
 Dim scale As Double
 If isScale Then
   scale = 1.0 / ( NUMBER + 1 )
 Else
   scale = 1.0
 EndIf

 ' For each of the sample outputs...
 Dim outerIndex As Integer
 For outerIndex = 0 To NUMBER
 
   ' Start sum at zero.
   outputData( outerIndex, 0 ) = 0.0
   outputData( outerIndex, 1 ) = 0.0

   Dim innerIndex As Integer

   ' For each input sample...
   For innerIndex = 0 To NUMBER
   
     ' Angle (radians) at this index.
     Dim angle As Double
     angle = multiplier * outerIndex * innerIndex / ( NUMBER + 1 )
     
     ' Sine and cosine of angle.
     Dim cosTerm As Double
     Dim sinTerm As Double
     cosTerm = cos( angle )
     sinTerm = sin( angle )
     
     ' Real and imaginary terms at this angle.
     Dim realTerm As Double
     Dim imaginaryTerm As Double
     
     realTerm = inputData( innerIndex + OFFSET, 1 ) * cosTerm
     realTerm = realTerm - inputData( innerIndex + OFFSET, 2 ) * sinTerm
     
     imaginaryTerm = inputData( innerIndex + OFFSET, 1 ) * sinTerm
     imaginaryTerm = imaginaryTerm + inputData( innerIndex + OFFSET, 2 ) * cosTerm
     
     ' Accumulate results.  
     outputData( outerIndex, 0 ) = outputData( outerIndex, 0 ) + realTerm
     outputData( outerIndex, 1 ) = outputData( outerIndex, 1 ) + imaginaryTerm
   Next
   
   ' Scale results.
   outputData( outerIndex, 0 ) = outputData( outerIndex, 0 ) * scale
   outputData( outerIndex, 1 ) = outputData( outerIndex, 1 ) * scale
 Next
 
 dft_guts = outputData
 
End Function

' Discrete Fourier Transform (DFT) in forward direciton (i.e. time->frequency).
Function dft( inputData As Variant )

 dft = dft_guts( inputData, -1, false )

End Function

' Inverse Discrete Fourier Transform (iDFT) in reverse direciton (i.e. frequency->time).
Function idft( inputData As Variant )

 idft = dft_guts( inputData, 1, true )

End Function
   To put this into a Calc document, goto Tool->Macros->Organize Macros->LibreOffice Basic.  Then you can assign the macro to a specific document, or put it in My Macros.  You will probably have to adjust security permissions just to allow the macro to run by going to Tools->Options->Security->Macro Security and setting the level from High (default) to Medium.  If you are not an idiot that shouldn't be a problem.
   This implementation isn't fast, but works fine for small data sets of around 50 samples.  Larger sets will still work, but you may be waiting awhile.

December 10, 2016

Implementing a Discrete Fourier Transform

Sitting fireside at my local coffee shop with a heavy snowstorm taking place outside, listening to Jehan Ariste Alain's Litanies, I came much closer to my goal of understanding how a Fast Fourier Transforms (FFT) works by understanding how to implement a Discrete Fourier Transform (DFT). The equation for a Fourier Transform:

Where g( t ) is some function with respect to time, f is some frequency, i is the square root of negative one. You can basically ask the function, what is the phase and amplitude for a given frequency for the time domain function g.

This is fine for transforming functions, but for most programming application we need to transform based on discrete data points. For this we can use a Discrete Fourier Transform (DFT):

Where xn is a data sample at time n, N is the total number of data samples, f is some integer frequency between 0 and N – 1, and G( f ) is amplitude of the given frequency. The discrete form is similar to the general form, except that the function g( t ) is replaced by our discrete sample, and there is now a finite range.

I have seen this equation in many places, but what the exponent is doing isn’t terribly obvious until you remember Euler’s formula:

Substitute this back in and we get:

Now we just have the imaginary number i as part of the equation, but this is fine because what we actually have is a complex number. They are generally represented in the form:

And here that breaks down pretty straight-forward:

In the form I’ve written it, G is a function. But it is actually a set just like x. So a slightly better form maybe:

There data, x, is actually complex data. Since the discrete data we are putting in consists of real values, it means the imaginary part is just 0. Now we just need to remember how to multiply complex numbers.

Where a and b are complex numbers made for a real part (nr) and an imaginary part (ni) represented n=(nr , ni).

With that, I actually have enough information to implement a Discrete Fourier Transform.

#include <math.h>
#include <stdio.h>

#ifndef M_PI
#define M_PI 3.1415926535897932384626433832795028841971693993751058209749445923
#endif

int main()
{
  enum { N = 16 };
  double inReal[ N ];
  double inImag[ N ];

  double outReal[ N ];
  double outImag[ N ];

  // Generate a signal.
  for ( unsigned index = 0; index < N; index += 1 )
  {
    inReal[ index ]  =       cos( 7 * 2 * M_PI * index / N );
    inReal[ index ] += 0.3 * sin( 3 * 2 * M_PI * index / N );
    inReal[ index ] += 0.5 * sin( 2 * 2 * M_PI * index / N );
    inImag[ index ] = 0;
  }

  // For each of the sample outputs...
  for ( unsigned k = 0; k < N; k += 1 )
  {
    // Start sum at zero.
    outReal[ k ] = 0;
    outImag[ k ] = 0;

    // For each input sample...
    for ( unsigned n = 0; n < N; n += 1 )
    {
      double cosTerm = cos( -2 * M_PI * k * n / N );
      double sinTerm = sin( -2 * M_PI * k * n / N );

      outReal[ k ] += inReal[ n ] * cosTerm - inImag[ n ] * sinTerm;
      outImag[ k ] += inReal[ n ] * sinTerm + inImag[ n ] * cosTerm;
    }
  }

  // Print the results.
  for ( unsigned index = 0; index < N; index += 1 )
    printf( "%u: %f %fi\n", index, outReal[ index ] / N, outImag[ index ] / N );

  return 0;
}

Here we construct a waveform consisting of 3 sine waves. One at 7 Hz, amplitude of 1, and phase 90°; one at 3 Hz, 0.3 amplitude, and 0°; and 2 Hz, 0.5 amplitude, and 0°. Then the transform, which consists of two loops. The outer loop is for each available integer frequency. The inner loop is the summation above. The last part prints the normalized of the transform, which agree with what I see when I do an FFT on the data.

This is a major step forward in my goal to understand an FFT. I still have a couple of questions, like why the output seems to be split between the upper and lower half. But I have proven (to myself anyway) that I can implement a DFT and that is a lot of progress.