Andrew Que Sites list Photos
Projects Contact
Main

January 12, 2017

Orthogonal Function with Dot Product Test

Playing some more with orthogonal functions to try and get a better handle of why a Fourier transform works. I decided I would try some simple functions I can work with easily to test orthogonality. The definition for as orthogonal function is one such that:

So I decided to use two very simple functions over the interval b = -a.

Orthogonal test:

Remember that ab even power b is the same for a and -a. Thus:

So my two functions are orthogonal. Now for something crazy.

I want to extract the coefficient for b. So I can to a projection for x2.

Notice how c does not show up in this result.

The same for x3.

Here, b does not show up in the result.

I verified this using a discrete check in a speadsheet. While I expected it would work, watching it actually do so was interesting.

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 29, 2016

Understanding the heart of a Discrete Fourier Transform

I have been working on understanding the innards of the Discrete Fourier Transform, and I’ve stumbled on what I believe to be the heart of the algorithm. Since then I have been trying to understand why it works. Let’s take a look at this heart:

These two functions will extract the amplitude and phase of the selected frequency c from the function f. With some experimentation, I have found the following to be true:

What this means is that if we have a sine wave at frequency b, amplitude a and phase d, we can extract phase and amplitude for some test frequency c. Anytime that b and c are equal, we get a non-zero value relating to phase and amplitude. If they are not equal, the answer is simply zero. This is precisely what is needed in a discrete Fourier transform, and why I call it the heart. In a DFT, the function checks the data for each possible integer frequency for the data set, and gets either it’s amplitude and phase, or zero.

First, the discrete Fourier transform in standard notation:

 

Use Euler's formula:

Now, assume x is real and break this into the real and imaginary parts and all frequencies:

 

Let’s try this on a set of data.

Here we have a (b) of frequency of 5 Hz with an amplitude (a) of 1 and phase angle (d) of 45°, and a check frequency (c) of 5 Hz. The sum of the real and imaginary parts is non-zero.

Now let’s try a check frequency of 6 Hz.

Here, the sum of the check frequencies for both real and imaginary are zero. It is more obvious when the check frequency is higher, say 50:

Here it is more obvious both the real and imaginary parts spend the same amount of time above and below 0, so overall summing to 0. It also doesn’t matter how many frequencies are in the signal. In the above examples there is just a single frequency, but it will extract amplitude and phase for the check frequency regardless of how many other frequencies are present in the signal.

I haven’t yet figured out why this is the case, but it makes perfect sense why this makes a DFT 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.