我没有声称这是最佳的，我也不认为你提到的过冲一定是不正确的。没有任何东西可以瞬间改变，因此当您从一个正弦曲线过渡到下一个正弦曲线时，您一定会看到一定范围的频率。

我对最小二乘估计的第一个想法是在展开阶段使用迭代器。由于频率是相位的导数，如果我们将频率积分，我们将得到相位。就线性系统而言，这看起来像

$$\mathbf{Ax} = \boldsymbol{\phi},$$

在哪里

$$\mathbf{A} = \begin{bmatrix} 1 & 0 & 0 & \cdots \\ 1 & 1 & 0 & \cdots \\ 1 & 1 & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}$$

是展开阶段。要获得展开的相位，您需要通过 I/Q 解调过程对实值信号进行分析。的最小二乘解将是频率的估计。$\boldsymbol{\phi}$$\mathbf{x}$

我拼凑的八度音阶代码如下。由于问题很小，我明确地形成了，但对于更大的问题，您可能希望使用类似的函数句柄来解决系统问题。我还在频率上运行了一个移动平均线，这样从一个到另一个的跳跃会更平滑一些，这应该会减少过冲。$\mathbf{A}$lsqr()

randn( 'seed', 71924 );
rand( 'seed', 7192 );

fmin = 100;
fmax = 1e3;
fs = 1.2 * 2 * fmax;
T = 5;
t = ( 0 : 1/fs : T ).';
Ns = length(t);

% Define some noise.
snr = 10.^( 50 / 10 );
n = sqrt( 1/snr ) * randn( Ns, 1 );

% Choose some random frequecies.
Nf = 50;
f = fmin + ( fmax - fmin ) * rand(Nf,1);

% Choose the intervals.
tmp = randperm( Ns );
inds = sort( tmp( 1 : Nf - 1 ) );

% Define the frequency for each sample.
fvec = zeros( Ns, 1 );
last = 1;
for ii = 1 : length(inds)
  fvec( last : inds(ii) ) = f(ii);
  last = inds(ii) + 1;
end
fvec(last+1:end) = f(end);

% Moving average over the frequencies to prevent
% bandwidth when the frequency changes.
fvec = conv( fvec, 1/3 * ones(3,1), 'same' );

% Define the sinusoid with time-varying frequency.
% We do so in a way that makes it phase continuous.
% We also add the noise here.
x = cos( 2*pi * cumsum( fvec ) ./ fs ) + n;

% The frequency sample locations.
df = fs / Ns;
floc = ( 0 : df : ( fs - df ) ) - ( fs - mod( Ns, 2 ) * df ) / 2;

% Shift the center of the positive half of the
% spectrum to zero.
y = x .* exp( -1j * 2*pi * t * fs/4 );

% Design a half-band filter.
Ntap = 51;
N = Ntap - 1;
p = ( -N/2 : N/2 ).';
s = sin( p * pi/2 ) ./ ( p * pi + eps );
s( N/2 + 1 )= 1/2;          
win = kaiser( Ntap, 6 );
h = s .* win;

% Form an analytic signal.
z = conv( y, h, 'same' ) .* exp( 1j * 2*pi * t * fs/4 );

% The unwrapped angle.
phi = unwrap( angle( z ) );

% Use an intergrator to get a least-squares
% estimate of the phase.
A = tril( ones( Ns, Ns ) );
fest = ( A \ phi ) * fs / (2*pi);

figure();
hold on;
plot( t, fest );
plot( t, fvec );
legend( 'Estimated', 'Truth' );
xlabel( 'Time (s)' );
ylabel( 'Frequency (Hz)' );