1
\$\begingroup\$

I'm trying to understand how to take the results of an FFT calculation and bin them to one-third octave bands. As my questions aren't about the code but about fundamentally understanding the arguments to and results from FFT, I assume this question fits better in electronics than in stackoverflow.

I start by calculating the center point and lower/upper range of each one-third octave band in Python:

freqs = [10**(0.1 * i) for i in range(12, 43)] # center frequency for 31 bands
fd = 10 ** 0.05 # 1.12201845
freqs_l = [freq / fd for freq in freqs] # lower
freqs_u = [freq * fd for freq in freqs] # upper

The minimum (freqs_l[0] == 14.125375446227551) and maximum (freqs_u[-1] == 17782.794100389234) values here give me the overall range of signals I'm interested in. Then I read in samples from my 16-bit, 48kHz WAV file using hound and run rustfft's forward-FFT on one second (48000 samples) of data:

const SAMPLE_RATE : usize = 48_000;
const NUM_BANDS : usize = 31; // this is the len(freqs) of the Python freqs above

let mut planner = FFTplanner::new(false);
let fft = planner.plan_fft(SAMPLE_RATE);

let reader = hound::WavReader::open(filename).unwrap();

let mut signal = reader
    .samples::<i16>()
    .take(SAMPLE_RATE) // just take one second of data for now
    .map(|samp| num::complex::Complex::new(samp.unwrap() as f32, 0f32))
    .collect::<Vec<_>>();

assert_eq!(signal.len(), SAMPLE_RATE);

// Convert time-domain signal to frequency-domain spectrum
let mut spectrum = signal.clone();
self.fft.process(&mut signal[..], &mut spectrum);

My main question here is around the frequency granuarlity of spectrum and how I can bucket these values into one-third octave bands. Both signal and spectrum are of size 48000, but does this cover the full Nyquist frequency range (0 - 24000 Hz), making each value 0.5Hz in width? For each input, I set an imaginary value of 0 (because all examples I've seen do this), but I have no idea if this affects the output in any important way.

Trying to find info on this, I read this SO question on finding the average FFT value of each frequency band and saw some code that pointed to numpy.fft.fftfreq that generates FFT frequencies for the sample length and sample spacing. From the former, It looks like I would calculate:

calculate d.

If I plug this into NumPy, numpy.fft.fftfreq(n=48000, d=0.14129655308810293), I get 48000 values, the first half non-negative (because it's symmetrical); the max value here is 3.538; I was expecting it would be closer to the Nyquist frequency.

I'm sure I'm way off here. How should I understand and adjust the FFT output frequencies?

\$\endgroup\$
3
  • \$\begingroup\$ No, "sample spacing" is just inverse sample rate. You've told fftfreq that you're taking a length-48000 FFT on a signal that's sampled at about 7.077Hz, and it's told you that the bins would then run from 1.4 millihertz to 3.538 hertz, but that's garbage in garbage out. \$\endgroup\$
    – hobbs
    Nov 23, 2020 at 7:28
  • \$\begingroup\$ Sanity check. If there's any doubt about interpretation of the theory, or about the software, FFT a known input frequency and see which bins contain energy. \$\endgroup\$ Nov 23, 2020 at 12:47
  • 1
    \$\begingroup\$ This doesn't address your main question, but you might want to look at the numpy function rfft which is slightly more efficient for real-valued time-domain signals (it doesn't require you to save the imaginary pats of a bunch of complex numbers whose real parts all are zero, and it only outputs the positive part of the spectrum because the negative part is known to be just the c.conjugate of the positive part) \$\endgroup\$
    – The Photon
    Nov 23, 2020 at 17:14

2 Answers 2

1
\$\begingroup\$

For an FFT, your number of data points should be a power of 2. You don't need to capture for a full second, in fact, your update speed will be really slow if you capture for 1 second. If this is a real-time display, the values below give an update speed of 3 Hz, still very slow.

You need to make some compromises. If this is a fun spectrum analyzer, it is not nearly as critical as if you were making a scientific instrument. The lowest bands will have some quantization error. If you are displaying in decibels, this won't be as much as you might think.

My fun spectrum analyzer only performs a 512 point FFT, I had to make a lot of compromises.

enter image description here

enter image description here

\$\endgroup\$
1
0
\$\begingroup\$

:-) to start from the end, in your formula you have calculated the number of decades (decimal orders) needed to span... not your desired audio spectrum, but instead the theoretical bandwidth of the digital recording, which means that you'd waste some of the bins on subsonic frequencies that no audio gear can reproduce anyway... and you divided that by the desired number of "bins of logarithmic frequency scale on your screen". So if you know that you want 31 log frequency bins on display, I'd calculate these from 20 kHz / 20 Hz = insert about 1000 instead of 24000 into the log10(). That's clearly 3 decades, /30 is about 10 bins per decade. Note 1000x is also about 10 octaves (2^10 = 1024) so 30 is indeed a pretty good number of log.bins, if you're aiming for 3 bins per octave.

Now... the FFT optimally takes 2^n time-domain samples and produces 2^n frequency lines (bins). If you have real-only input (no imaginary component = not complex), I believe the output will only contain meaningful data in the lower half of the bins (which will still be complex numbers), and the "bins above Nyquist" should optimally contain zeroes - as a proof that you have filtered your analog input enough, which you probably won't be able achieve 100% (i.e. under the level of your ADC's quantisation noise, to steer completely free of mirror images).

The complex output bins mean "amplitude and phase shift". Before you try to take an average, I suggest that you first calculate an absolute value of each complex output bin, to get the amplitudes (and ignore the phase information). After that, note that the FFT produces evenly spaced output bins (2^n of them), and you want a logarithmic scale for display. So each of your log display bins will span a different number of raw FFT output bins - I'll leave it up to you if your average will contain some sort of logarithmic weighting, maybe even including fractions of "bins on the edge" (as the edges won't fit on integer count boundaries of your raw FFT output bins).

While you're playing with the maths, don't be afraid to be creative, but try to understand what you're doing :-) I wish you good luck.

\$\endgroup\$
2
  • \$\begingroup\$ I read your response several times and played with the code; it sounds like I need to get my frequency bins as bins = numpy.fft.fftfreq(1000, 1/48000). This gives me 1000 bins with 48Hz between each. Using the first half (which are non-negative), this is a 24000 range. The last one covers ~76.2 × 48Hz bands from (14125, 17782), so I need to add the 76 FFT output values to cover that one bin. Since I only need freqs from 1.58Hz to 15.8kHz, I then don't need all 48000 FFT output values - and these first half are all I can use anyway if I'm using the first 500 fftfreq bins. Is that right? \$\endgroup\$
    – user655321
    Nov 24, 2020 at 17:17
  • 1
    \$\begingroup\$ @user655321: I didn't re-do all the maths, but it sounds about right, yes :-) Mattman944 has uploaded a neat spreadsheet in his response. \$\endgroup\$
    – frr
    Nov 25, 2020 at 14:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.