如果我们以足够的采样率周期性地对连续带限信号进行采样，我们可以通过使用离散傅里叶变换 (DFT) 来估计其频谱。有限离散信号样本的 DFT 可找到一组预定义频率的幅度：每个频率都是基频的倍数。
但是假设原始连续信号主要由不可通约的频率组成，例如。然后，即使只有少数这样的频率，我们样本的 DFT 实际上在每个点上都将是非零的，因为组成频率将错过 DFT 所代表的频率，并且它们的峰值会变宽。
现在，为了获得更高分辨率的频谱估计，我们可以增加采样时间，从而使基频变小。但是要对频率进行非常好的估计，这可能意味着采样时间增加一百倍或更多，这并不总是可行的。
似乎对于不太丰富的频谱，我们应该能够从原始的小样本中找到更好的频率集。例如，对于单频信号，我们可以尝试最小化来自测试波的信号差异能量，在频率和相位的（连续）空间中寻找最小值。类似的过程可以应用于少数频率的信号，在频率上的每次迭代中从中移除频率分量的最佳估计。
上述最小化过程似乎是一种蛮力解决方案。有没有更有效的方法来实现同样的目标？
查找不相关频率的常用方法是如您所建议的：找到最大的正弦曲线，将其从信号中删除，然后找到下一个最大的，将其删除，等等。
下面是一些执行此操作的 matlab 代码。
首先，生成一个三正弦噪声信号。这是在这里绘制的：
然后，我只使用Quinn-Fernandez 估计器估计最大正弦波的频率（它既便宜又令人愉快），每次都去除正弦波。
这产生了下面的三个图。注意最后一个基本上只是噪音。
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)));
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
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
%
% Copyright (C) 1993 CRC for Robust & Adaptive Systems
%
% 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.