Accurate frequency estimation with closely-spaced tones

Question

I have a MIPs-limited DSP platform where I am taking an FFT of a streaming audio signal at regular intervals. I need to find the top 4 amplitude peaks and their corresponding frequency for every frame. I am using the "standard" FFT interpolation formula to get accurate frequency estimation when the true peak is between 2 bins, and this normally works very well. Because I am so MIPs-constrained, I can only do a 1K fft at fs=48KHz. My problem is that when the input contains low-pitched male speech, the fundamental tone is around 100 Hz (bin 2), and the 2nd-harmonic is around 200hz (bin 4), so there is only 1 bin that lies between the two peaks. That means that the "shoulder" bins used in the interpolation formula contain leakage energy from both both components, and this reduces the accuracy of the interpolation by quite a large factor. I know the "right" answer is to use a longer FFT, but that is just not an option in this case. Any ideas? The solution needs to be simple, or it likely won't fit on my DSP.

Answer 1

Here's a GNU Octave (Matlab clone) script that demonstrates using the cross spectrum to turn individual FFT bins into FM detectors to get estimates of the instantaneous frequency in each bin:

Fs = 48000; % samples/sec
N = 512; % FFT size

% Test signal is a sum of cosinusoids
f = [110, 220].'; % Hz.  Test signal frequencies
phi0 = 2*pi*rand(length(f),1); % radians.  Test signal initial phases

M = 8; % Number of FFT blocks in test signal
L = M * N + 1; % Number of samples in test signal
k = [0:(L-1)]; % sample index

% Build the test signal
x = ones(1,length(f))*cos(2*pi/Fs*diag(f)*repmat(k,length(f),1)+phi0);

% Compute the cross-spectra for the FFT blocks we have in the test signal
pkg load signal % Octave requires for Blackman-Harris window
%W = blackmanharris(N).';
%W = kaiser(N,3.0).';
W = rectwin(N).';
cross_spectra = zeros(M, N);
for j = [0:(M-1)]
    n0 = j*N + [1:N];
    X0 = fftshift(fft(x(n0).*W,N));
    X1 = fftshift(fft(x(n0+1).*W,N));
    cross_spectra(j+1,:) = X1.*conj(X0);
end

% Argument of the cross spectrum approximates dphi_dt for each bin
bin_dphi_dt = arg(cross_spectra);

% Which is scaled to an approximate instantaneous frequency for
% an assumed single frequency signal component in that bin
inst_freq_est = bin_dphi_dt * Fs/2/pi;

% A simple filter to smooth out the estimate.
avg_inst_freq_est = mean(inst_freq_est);

bins = [(-N/2):(N/2-1)];
freqs = bins/N*Fs;

% Get the frequency bins of the known freqs in our test signal
[f_idx, bin_idx] = find((f > (freqs - Fs/N/2)) .* (f < (freqs + Fs/N/2)));

printf("Freq values: Bin Center: %f, Truth: %f , Estimated: %f\n", [freqs(bin_idx).', f, avg_inst_freq_est(bin_idx).'].')

figure(1);
plot(freqs, 10*log10(mean(abs(cross_spectra))), '-x');
title('Average Cross Spectrum Power');
xlabel('Frequency (Hz)');
ylabel('Cross Power (dB)');

figure(2);
plot(freqs, avg_inst_freq_est, 'x');
title('Cross Spectrum Phase: Spectrum Phase Velocity Average scaled to Instantenous Frequency');
xlabel('Frequency Bin Center (Hz)');
ylabel('Instantaneous Frequency Estimate Average (Hz)');

A few sample runs of the program yield this output:

>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 108.306172
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 215.834249
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 112.496531
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 214.549783
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 106.189949
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 210.683539
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 106.526426
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 212.488463
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 114.102267
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 226.823594
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 113.526806
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 223.008182
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 111.629428
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 221.790921
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 111.100902
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 217.472948
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 108.005589
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 222.292734
>> cross_power_freq_est
Freq values: Bin Center: 93.750000, Truth: 110.000000 , Estimated: 116.300657
Freq values: Bin Center: 187.500000, Truth: 220.000000 , Estimated: 232.532066