FFT (Fast Fourier Transform)

FFT is a mathematical tool used to decompose a signal to its underlying base frequencies (known as harmonics). This is very helpful to detect the periodicity aspect in a signal.

To use FFT in greycat, you need to add the algebra library and use transforms module as following:

@library("algebra");
use transforms;

FFT class

In order to show a usage example of FFT, let’s first generate a timeseries at random and store it in a nodeTime as following:

fn main() {
    var ts = nodeTime::new();
    //Random number generator with a fixed seed
    var rand = Random::new();
    rand.setSeed(1234);

    var starting = time::now().to(DurationUnit::microseconds);
    var avgT = 10000; //generate a timepoint every 10 ms
    var N = 17;

    var offset = 0;
    println("Displaying original timeseries");
    for (var i = 0; i < N; i++) {
        offset = offset + avgT;
        var v = rand.uniformf(0.0,100.0);
        ts.setAt(time::new(starting + offset, DurationUnit::microseconds), v);
        println("${i} - ${offset} => ${v}");
    }
    println("");
}

Next step, let’s store some basic information about this timeseries in variables, add the following code in the end of the previous main:

    var ts_start = ts.firstTime()!!; //Get the first timepoint of the timeseries
    var ts_end = ts.lastTime()!!;   //Get the last timepoint of the timeseries
    var size = ts.rangeSize(ts_start, ts_end); //get the total number of elements in the timeseries

Now we have the start time, end time and the size of the timeseries. FFT algorithm is optimized to run better on specific sizes, so to get the optimize size run the following:

    var best_size = FFT::get_next_fast_size(size); //gets an optimized size for FFT
    var sampling_step = (ts_end - ts_start).to(DurationUnit::microseconds) / (best_size - 1);
    var sampling_duration = duration::new(sampling_step, DurationUnit::microseconds);
    info("time series has ${size} timepoints => next best size is: ${best_size}, sampling duration is: ${sampling_step} us");

To create a forward transform FFT, we need to create the following:

    var fft = FFT::new(best_size,false); //Creates a forward transform fourier
    var time_complex = Tensor {};   //Create a time complex tensor to store the timeseries

Note that the FFT object is un-serializable, so it can’t be stored to disk. It contains however an array of size N of precalculated cosines (to accelerate the computation). So if you need to run fourier transform over windows of same length, think to re-use this object for performance.

Running forward transform

To fill the time_complex tensor, we need to iterate on the nodeTime at a fixed sampling rate, we do so using the following:

    //Initialize the time_complex with c128 tensor type, and dimension best_size
    time_complex.init(TensorType::c128, [best_size]); 

    //Initialize a nodeTimeCursor to iterate on the nodeTime
    var nc: nodeTimeCursor = nodeTimeCursor::new();
    nc.init(ts);

    println("Displaying sampled timeseries");
    for (var i = 0; i < best_size; i++) {
        //set the time_complex tensor 
        time_complex.set([i], nc.current() as float);

        println("${i} - ${nc.currentTime()!!.to(DurationUnit::microseconds) - starting} => ${nc.current()}");
        nc.skip_duration(sampling_duration);
    }
    println("");

Last step, in order to convert time_complex tensor to frequency_complex we do the following:

    //Convert time complex to frequency complex
    var frequency_complex = Tensor {};
    fft.transform(time_complex, frequency_complex);

After the fft.transform call you get the frequency tensor as a complex object containing the FFT decomposition of the time_complex tensor.

Frequency table

FFT class contains several helpers to improve the usage of the fft_complex tensor in gcl. One such helper is the get_frequency_table:

    var frequency_table = FFT::get_frequency_table(frequency_complex,sampling_duration);

This method converts the frequency_complex tensor into a table useful for display in the frontend. The table contains the following columns:

Low pass filter

    var lowpass = FFT::get_low_pass_filter_size(frequency_complex, 0.95);
    println(lowpass);

This method gives the index of where to cut the frequency spectrum to keep 95% of the signal. In signal processing it’s called low pass filtering.

FFT model

If you don’t want to pilot the FFT by yourself, and you care about getting a quick FFT model of a timeseries, we have the easy to use following methods:

var fft_model = FFTModel::train(ts, ts.firstTime()!!, ts.lastTime()!!);

This creates automatically an FFT decomposition of a timeseries, you just need to provide the timeseries, first time-point where to start the analysis and last time-point. Inside the model you will find the previously discribed best_size, frequency_complex, frequency_table etc

Extrapolate method

This fft_model can be used to generate a table of timeseries from any time-point to any time-point as following:

var result = fft_model.extrapolate(start_time, end_time, low_pass_filter_ratio, nb_harmonics, skip_elements);

the low_pass_filter_ratio, nb_harmonics, skip_elements are not mandatory attributes.