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.
-
$\begingroup$ You can probably use the cross spectrum - multiplication of the current FFT with the complex conjugate of a previous time FFT - to home in on the peak. The cross spectrum turns each freq bin into an FM detector really. Look for papers by Douglas Nelson. $\endgroup$– Andy WallsDec 22, 2022 at 0:40
-
$\begingroup$ Thanks Andy, I’ll look that up! $\endgroup$– BobDec 22, 2022 at 6:16
1 Answer
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
-
$\begingroup$ This is great, thanks for the detailed response ! $\endgroup$– BobDec 24, 2022 at 13:20
-
$\begingroup$ @Bob You're welcome. I hope it works for you. Since it's 2 FFTs per set of raw frequency estimates, I knocked the FFT size down to 512 bins, and hence the bin size is about 94 Hz wide. So if you are unlucky enough to have two frequencies closer together than that in the same bin, you won't be able to discriminate them as 2 tones. Whereas you may have at least known you had 2 tones like that with the 1024 point FFT, even if you couldn't have estimated their frequencies well. $\endgroup$ Dec 24, 2022 at 21:44
-
$\begingroup$ BTW, the fftshift() operations in the Octave script are totally unnecessary. I just used them for plotting graphs. $\endgroup$ Dec 24, 2022 at 21:46