Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Some example code #48

Open
larsmob opened this issue Mar 31, 2024 · 1 comment
Open

Some example code #48

larsmob opened this issue Mar 31, 2024 · 1 comment

Comments

@larsmob
Copy link

larsmob commented Mar 31, 2024

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
Copy link
Owner

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants