liamappelbe / fftea

A simple and efficient FFT implementation for Dart
Apache License 2.0
63 stars 1 forks source link

Some example code #48

Open larsmob opened 5 months ago

larsmob commented 5 months ago

I wrote a small bit of code to make sure that I understand how this package works and that I get correct results. Maybe this is useful for someone else.

import 'dart:math';
import 'package:fftea/fftea.dart';
import 'dart:typed_data';

const double time = 1; // Seconds
const int fs = 48000; // Samplerate
const double delta_t = 1/fs; // Time between each audio sample
const double delta_f = fs/(time*fs); // Frequency resolution of fft output
const double ampl1 = 0.8, freq1 = 256, angle1 = 45; // Sine wave parameters (for generating input waveform)
double log10(num x) => log(x) / ln10;

void myfft(List<double> tvec, [List<double>? window, double windowcorrectionfactor = 1]) {
  List<double> itvec = List.from(tvec);
  if (window != null) {
    // Apply window
    for (var i = 0; i < itvec.length; i++) {
      itvec[i] *= window[i];
    }
  }
  // Calculate fft
  final fft = FFT(itvec.length);
  final freqbins = fft.realFft(itvec);
  // Remove frequencies above the Nyquist frequency
  Float64x2List fftresultF64x2 = freqbins.discardConjugates();
  // Calculate the phase/angle of each frequency component.
  List<double> fftangle = fftresultF64x2.map((i) => atan(i.y/i.x).toDouble()).toList();
  // Get the magnitude of each frequency component.
  Float64List fftamplF64 = fftresultF64x2.magnitudes();
  // Normalise (multiply by two because we discarded half the energy, divide by number of input samples) and correct for Window function
  List<double> fftresult = fftamplF64.map((i) => windowcorrectionfactor*2*i.toDouble()/itvec.length ).toList();

  // Find max, this should be from our input sine wave
  double max = 0;
  int maxindex = 0;
  for (int i = 0; i < fftresult.length; i++) {
    if (fftresult[i] > max) {
      max = fftresult[i];
      maxindex = i;
    }
  }
  // Look at bins around the max bin.
  for (int i = maxindex-3; i <= maxindex+3; i++) {
    //print("bin: ${i}, value: ${fftresult[i]}, ${freqbins[i].x}, ${freqbins[i].y}, ${freqbins[48000-i].x}, ${freqbins[48000-i].y}");
  }

  print("Number of input samples: ${itvec.length}, output bins: ${freqbins.length}, after discard: ${fftamplF64.length}");
  print("Amplitude is ${max} and phase is ${360*fftangle[maxindex]/(2*pi)} at frequency ${maxindex*delta_f}");
  print("");

}

void fftstuff() {
  List<double> tvec = []; // A list for holding the input samples
  // Generate sine to feed the fft.
  for (double i = 0; i <= time; i += delta_t) {
    tvec.add(ampl1 * cos(2*pi*freq1*i+angle1*2*pi/360));
  }

  // Perform fft with no window (rectangular), then hanning and hamming.
  myfft(tvec);
  myfft(tvec, Window.hanning(tvec.length), 2);
  myfft(tvec, Window.hamming(tvec.length), 1.85);
}

void main() {
  fftstuff();
}

The output is:

Number of input samples: 48000, output bins: 48000, after discard: 24001
Amplitude is 0.800000000000011 and phase is 45.000000029805 at frequency 256.0

Number of input samples: 48000, output bins: 48000, after discard: 24001
Amplitude is 0.7999833333354626 and phase is 45.00000002217547 at frequency 256.0

Number of input samples: 48000, output bins: 48000, after discard: 24001
Amplitude is 0.7991858166684803 and phase is 45.000000023305795 at frequency 256.0
liamappelbe commented 1 month ago

Sorry for the slow reply. This repo could use more examples, so feel free to clean this up and send a PR. One note I have is that windowing is mainly useful for STFT. It's not usually used for a one shot FFT like this.