这里有一些各种频率估计器的 Matlab 实现。 但是，它们可能会被高次谐波“欺骗”。

你知道一个范围，也许最好的技术是 Eric Jacobsen 的估计器之一。看看他的页面。

本质上，Eric 的估计器在 X(k) 处使用 FFT bin，在 ) 的任一侧为峰值。如果您选择 FFT 长度以匹配您的已知范围（以便和和，那么你应该很高兴。 $X(k)$$X(k+1), X(k-1)$$k \rightarrow 3.45$$k+1 \rightarrow 3.50$$k-1 \rightarrow 3.4$

谐波乘积频谱 (HPS) 似乎对您有用，HPS 算法的一个特点是将多个谐波相乘，这会为您提供基频的峰值，我相信您需要像 12288 这样的大 FFT 点来采样信号在 44100hz，因为它的第一个 bin 是 3.5889hz。HPS 源代码示例： https ://stackoverflow.com/questions/19765486/matlab-code-for-harmonic-product-spectrum

如果谐波恰好是谐波（并且不是许多物理系统产生的轻微不谐波），并且您对强谐波有良好的频率估计，使得该谐波频率估计适合基波搜索范围的唯一整数倍频率，那么您可以只除以该整数倍以获得基频的频率估计。

一旦有了频率估计，您就可以使用 1 个 DFT 或 Goertzel 滤波器计算幅度估计，其长度恰好是该频率估计周期的整数倍。（...或 3 个非常接近的长度，然后是抛物线峰值插值。）

非常感谢，这看起来很有希望。

我使用 Erics 代码（Macleod 算法）来实现以下片段：

% Amplitude of signal.
ACW=1;

% Specify signal.
Fs = 8000;                   % samples per second
dt = 1/Fs;                   % seconds per sample
nHarmonic = 7;               % order of harmonic
Fc = 21.5;                   % Hz
StopTime = 0.25;             % seconds
t = (0:dt:StopTime-dt)';     % seconds

% Specify evaluation tolerance
Tol = 0.01;

% Effective length of test vector.
tstlen=length(t);

% Generate signal
cw = ACW * sin(2*pi*(1:tstlen)*Fc/tstlen) + 1.5 * ACW * sin(2*pi*(1:tstlen)*(Fc/tstlen) * 3) ...
  + 2.5 * ACW * sin(2*pi*(1:tstlen)*(Fc/tstlen) * 5)+ 0.3 * ACW * sin(2*pi*(1:tstlen)*(Fc/tstlen) * 7);

% Number of trials in test.
N=10000;           

% Desired bin number of tone relative to long test length.
bin=Fc; 

% Allocate vector for results.
macldest=zeros(size(1:N));

% Calculate desired target result for comparison.
target=bin;

fprintf('Peak is at %f.\n',bin);

for I = 1:N,        % Run trials.
% Calculate signal power.
sigp=sum(abs(cw(1:tstlen).^2));         

% DFT.
dftshrt(1:tstlen)=fft(cw(1:tstlen));

% DFT magnitude.
magshrt(1:tstlen)=abs(dftshrt);

% Find raw peak magnitude and location.
[rawmag,rawind]=max(magshrt);           

% Isolated 3 samples around peak.
pk3vect(1:3)=dftshrt(rawind-1:rawind+1);    

%Do Macleod's estimation.
ref = pk3vect(2);                               % Isolate phase reference.
R = real(pk3vect.*conj(ref));                   % Generate phase corrected coefficient vector.
gamma = (R(1)-R(3))/((2*R(2))+R(1)+R(3));       % Calculate offset.
delta = (sqrt(1 + 8*gamma*gamma)-1)/(4*gamma);  % Final estimate.
macldest=rawind-1+delta;        
end

% Calculate average result.
macldr=mean(macldest)

% Evaluate estimated frequency
Factor = round(macldr/Fc)

if (macldr > (Fc-Tol) & macldr<(Fc+Tol) )
  result=macldr
elseif ((macldr/Factor)>(Fc-Tol) & (macldr/Factor)<(Fc+Tol))
  result=macldr/Factor
end

但是现在估计基本面峰值最准确的方法是什么？