# Is there an efficient algorithm to precisely find the spectrum of a finite discrete signal consisting of incommensurable frequencies?

If we periodically sample a continuous band-limited signal at a sufficient sampling rate, we can estimate its spectrum by using Discrete Fourier transform (DFT). DFT of a finite discrete signal sample finds amplitudes of a predefined set of frequencies: each frequency is a multiple of the fundamental.

But suppose the original continuous signal consists mostly of incommensurable frequencies, e.g. $$a\exp(it)+b\exp(\sqrt2 it)+c\exp(\pi it)$$. Then, even if there's only a handful of such frequencies, the DFT of our sample will be nonzero virtually at every point, because the constituent frequencies will miss the ones represented by the DFT, and their peaks will be widened.

Now, to get a higher-resolution estimate of the spectrum, we could increase sampling time, so that the fundamental frequency gets smaller. But to get a really good estimate of the frequencies this might mean a hundredfold or more increase in sampling time, which is not always practical.

It seems that for a not too rich spectrum we should be able to find better set of frequencies from the original, smallish sample. E.g. for a single-frequency signal we could try minimizing energy of difference of the signal from a test wave, searching for the minimum in the (continuous) space of frequencies and phases. Similar process could be applied to a few-frequencies signal, removing the best estimate of a frequency component from it at each iteration over the frequencies.

The above minimization procedure seems like a brute-force solution. Are there more efficient ways to achieve the same?

The usual way to find unrelated frequencies is as you suggest: find the largest sinusoid, remove it from the signal, and then find the next largest, remove it, etc.

Below is some matlab code that does this.

First, a three-sinusoid noisy signal is generated. This is plotted here:

Then, I just estimate the frequency of the largest sinusoid using the Quinn-Fernandez estimator (it's cheap and cheerful), removing the sinusoid each time.

This yields the three plots below. Note the last one is basically just noise.

# Matlab code only below.

a = 1;
b = 2;
c = 3;
N = 1000;
fs = 10000;
omega1 = 2*pi*1234;
omega2 = exp(1)*1234;
omega3 = sqrt(2)*pi*1234;
rng(1234);
phi1 = rand(1)*2*pi;
phi2 = rand(1)*2*pi;
phi3 = rand(1)*2*pi;

t = [0:N-1]/fs;

x = a*cos(omega1*t+phi1) + b*cos(omega2*t + phi2) + c*cos(omega3*t+phi3) + 0.1*randn(1,N);
figure(1);
clf;
subplot(211);
plot(t,x);
subplot(212);
plot(abs(fft(x)));

[omegaA, phiA, AmpA, x2] = removeSinusoid(x, fs);
[omegaB, phiB, AmpB, x3] = removeSinusoid(x2, fs);
[omegaC, phiC, AmpC, x4] = removeSinusoid(x3, fs);

figure(2);
clf;
subplot(311);
plot(abs(fft(x2)));
subplot(312);
plot(abs(fft(x3)));
subplot(313);
plot(abs(fft(x4)));


# Remove sinusoid function

function [omegaHat, phihat, Amphat, xremoved] = removeSinusoid(x, fs)

N = length(x);
omegaHat = qnf(x');
DFT = sum(x.*exp(-1j*omegaHat*[0:N-1]));
phihat = angle(DFT);
Amphat = abs(DFT)/N*2;

xremoved = x - Amphat*cos(omegaHat*[0:N-1]+phihat);
omegaHat = fs*omegaHat;
end


# Quinn-Fernandes Implementation

function [est] = qnf(signal)
% QNF - The Quinn - Fernandes frequency estimator.
%
%          Inputs:  signal - T x N matrix where
%                            T = data length
%                            N = number of signals
%                            (i.e. N signals in columns).
%
%          Outputs: est - N Quinn-Fernandes frequency estimates.
%
% [1]  B.G. Quinn & J.M. Fernandes, "A fast technique for the estimation of frequency,"
%      Biometrika, Vol. 78(3), pp. 489--497, 1991.

% $$Id: qnf.m 1.1 2000/06/07 18:57:16 PeterK Exp PeterK$$

% File: qnf.m
%
%
% This software provided under the GNU General Public Licence, which
% is available from the website from which this file originated. For
% a copy, write to the Free Software Foundation, Inc., 675 Mass Ave,
% Cambridge, MA 02139, USA.

%
% Initializations
%
[t,ns]=size(signal);
xb=mean(signal);
signal=signal-ones(t,1)*xb;
t3 = t+1;
vrsn = version;
if vrsn(1)=='3'
y=fft([signal; zeros(signal)]);
else
y=fft([signal; zeros(size(signal))]);
end

z=y.*conj(y);
z=z(2:t3,:);

[m,j]=max(z(2:t-1,:));

j=j+1;

a=2*cos(pi*j/t);
y=y(1:2:2*t,:);

%
% Quinn-Fernandes method
%
b=[1];
nm=t-1;
for jjj=1:2

for q = 1:ns
c=[1;-a(q);1];
y(:,q) = filter(b,c,signal(:,q));
end

v = sum(signal(2:t,:).*y(1:nm,:))./sum((y(1:nm,:).*y(1:nm,:)));
a = a+2*v;
end

est=acos(a/2);

% Author: SJS 1992; Adapted from code within ttinpie.m (author PJK)
%
% Based on: P.J. Kootsookos, S.J. Searle and B.G. Quinn,
% "Frequency Estimation Algorithms," CRC for Robust and
% Adaptive Systems Internal Report, June 1993.

• Interesting that Quinn-Fernandes can estimate the frequency of the largest sinusoid even when there are other sinusoids in the signal. Do you know how close the sinusoids can be (in frequency and in amplitude) before QN becomes inaccurate?
– MBaz
Commented Nov 5, 2021 at 17:24
• @MBaz I don't have a good feeling for how that varies. I suspect for like-amplitude tones that they're going to have to be separated by at least an FFT bin with ($2\pi/N$).
– Peter K.
Commented Nov 5, 2021 at 19:21
• Thanks! I'll run some experiments to see how it behaves.
– MBaz
Commented Nov 5, 2021 at 19:45