我们如何准确估计这样一个信号的频率？

Question

我有一个信号，其表达式如下所示：

当信号中的A、B、ω、D、E都是未知参数时，用什么方法可以准确拟合这个信号呢？特别是，其中最关键的未知参数是 ω。

信号看起来像这样：

如何准确估计其中最关键的ω的值？

我尝试了MATLAB的lsqcurvefit和nlinfit函数，但精度差强人意。尝试使用FFT提取ω，发现结果也不是很准确。

Answer 1

1.-问题所附图中提供的信号清楚地表明 一个单音 + 1 个锯齿波。

这是提供的“屏幕截图”类型

只要音调有足够的功率来显示上面锯齿波谱 然后估计单个音调的搜索频率足以发现 被锯齿波谱包围的单音峰值。

如果所寻找的音调频率发生在锯齿波包络以下的功率下，那么您必须沿着锯齿波频谱包络寻找异常，否则，在没有音调的情况下，随着频率的增加，会出现平滑的

1/k

衰减。

2.- 您正在分析的信号并不完全是所提供的信号 表情

这是提供的信号表达式

即使提供的信号是您所拥有的全部，为了使用傅里叶变换，也必须假设有限的功率信号，而不是无限的能量信号。

改写；信号

x1(t)=a*t

、

x2(t)=a*|t|

没有沿

[-Inf Inf]

或

[0 Inf]

进行傅立叶变换，因为频谱变换的积分不会产生任何结果。

因此我们需要对信号进行时间窗口处理，因此要考虑进行频谱分析的信号是 1 个音调和 1 个锯齿的总和。

由于没有提及

tstart

tstop

或

t

时间参考，我假设有一个 :

f1=20                       % [Hz] tone frequency
f2=1                            % [Hz] sawtooth base frequency

if f1>=f2
    dt=1/(50*f1);       % time step [s] seconds
    else
       dt=1/(10*f2);
end
Fs=1/dt;                    % [Hz] sampling frequency

tstart=0                    % start time [s]
tstop=5                     % stop time [s]
t=[tstart:dt:tstop];        % time reference
L=numel(t);

A1=.5                       % tone amplitude            
A2=2                        % sawtooth peak2peak amplitude

如果

f2*tstop

不是整数，则应添加

min(f1,f2)

的分数或删除最后一个循环分数。

为了简单起见，我没有包括需要处理循环尾部的情况。

% 1 tone
n2=f2*tstop;           
x1=A1*sin(2*pi*f1*t);

% sawtooth
T2=1/f2;
x20=[0:dt:T2];        % sawtooth base cycle
x2=A2*repmat(x20(1:end-1),1,n2);
x2=x2-A2/2;
y1=x1(1:min(numel(x1),numel(x2)))+x2(1:min(numel(x1),numel(x2)));
t=t(1:min(numel(x1),numel(x2)));

figure(1);
plot(t,y1);
grid on;
title('Sawtooth + Tone');
xlabel('t[s]');

所以很明显，这个问题中提供的信号图像是 1 个音调和 1 个锯齿波。

3.- 提供的信号图像不是 2 个音调和 1 个锯齿

包含相同频率的 2 个音调和 1 个锯齿波的信号实际上如下所示：

A2=A1;                   % A2 is B in the question
y2=y1+A2*cos(2*pi*f2*t);

figure(2);
plot(t,y2);
grid on;
title('Sawtooth + 2 Tones same amplitude');
xlabel('t[s]');

A2=2*A1;                   
y2=y1+A2*cos(2*pi*f2*t);

figure(3);
plot(t,y2);
grid on;
title('Sawtooth + 2 Tones 1:2 amplitudes');
xlabel('t[s]');

很明显，问题中提供的信号图像是 1 个音调和 1 个锯齿。

随着时间信号的产生，时间不连续性可以得到解决，因此只有平滑的过渡。不难，但是费力，而且为了回答这个问题，这样的努力不会增加相关内容。

4.- 光谱

X1=fft(x1);
XP12=abs(X1/L);
XP1=XP12(1:floor(L/2)+1);
XP1(2:end-1)=2*XP1(2:end-1);
f = Fs*(0:floor(L/2))/L;

figure(4);
plot(f,XP1);
grid on;
title('X1(f)=TF(x1(t)) single tone : single side spectrum');
xlabel('f');

与单音相比，锯齿波具有丰富的频谱、噪声和干扰

Y1=fft(y1);
YP12=abs(Y1/L);
YP1=YP12(1:floor(L/2)+1);
YP1(2:end-1)=2*YP1(2:end-1);
f = Fs*(0:floor(L/2))/L;

figure(6);
ax6=gca
plot(ax6,f,YP1);
grid(ax6,'on');
title(ax6,'Y1(f)=TF(y1(t)) tone + sawtooth : single side spectrum');
xlabel(ax6,'f');
hold(ax6,'on')

5.- 锯齿信号基础知识参考：

我们可以回顾一下锯齿信号的基础知识：

https://mathworld.wolfram.com/FourierSeriesSawtoothWave.html

和

https://en.wikipedia.org/wiki/Sawtooth_wave

6.- 锯齿信号频谱包络

T2

是锯齿基循环。

如果锯齿波的时间定义为无 DC，则删除 DC 常数。

A2

是锯齿波幅度 WITH DC 正如我最初定义的那样。

零直流锯齿波，幅度为

A2/2

。

锯齿信号（非功率）频谱的包络是

1/(pi*k)

，但这不包括

Fs

和

L

：

k1=[0:1:numel(f)-1];

% sawtooth signal (NOT POWER) spectrum envelope
S=2*L/Fs*A2/pi*1./(k1);       % double the signal (NOT POWER) envelope because it's 1-side spectrum
plot(ax6,f,S,'Color','r','LineWidth',2)

评论

注意没有 DC，因为音调是

sin

并且我已经删除了所有锯齿 DC。

对于 2 个音调，由于两者处于相同频率，因此幅度谱将相似。

如果您想了解如何捕捉锯齿信号杂乱的弱音，请发布另一个问题。