Bucketing FFT into one-third octave bands

Question

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:

.

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?

Answer 1

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.

Answer 2

:-) 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.