提问人:BlueConda 提问时间:8/7/2023 更新时间:8/31/2023 访问量:75
我们如何准确估计这种信号的频率?
How can we accurately estimate the frequency of such a signal?
问:
我有一个信号,其表达式如下所示:
当信号中的A、B、ω、D、E都是未知参数时,我可以使用什么方法来准确拟合该信号?特别是,其中最关键的未知参数是ω。
信号如下所示:
如何准确估计其中最关键的ω值?
我尝试了 MATLAB 的 lsqcurvefit 和 nlinfit 函数,但准确性很难令人满意。尝试使用 FFT 提取 ω,发现结果也不是很准确。
答:
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
如果不是整数,则应添加 a 的分数或削减最后一个周期的分数。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 个锯齿
包含相同频率和 1 个锯齿的 2 个音调的信号实际上如下所示:
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
是锯齿基循环。
如果锯齿被时间定义为无直流电,则删除直流常数。
A2
是我最初定义的带直流的锯齿振幅。
零直流锯齿,具有振幅。A2/2
锯齿信号(非功率)频谱的包络是,但这不包括和:1/(pi*k)
Fs
L
k1=[0:1:努梅尔(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,因为音调是,我已经移除了任何锯齿 DC。sin
对于 2 个音调,由于两者都在同一频率上,因此幅度谱相似。
如果您想了解如何捕捉带有锯齿信号的微弱音调,请发布另一个问题。
按照 SE 成员 Jean Jacquelin 的想法,假设数据很密集,很容易得到第一个估计值。
在这种情况下,u 函数 y 的双数值积分,即 SSy,具有 SSy = -w^2 y + D/6 t^3 + E/2 t^2 + U t + Z 的性质。 因此,使用简单的线性方程,我们得到 w。 有了已知的 w,整个问题就变成了纯线性的,我们可以拟合其他参数。通过将其代入非线性拟合,可以改善结果。不过,人们必须小心。正弦-余弦拟合在引入相位时效果更好(请参阅此处的一些详细信息)。 最后,我们可以尝试使用 FT 来提高频率。首先,我删除偏移量和斜率。为了精确确定频率,我使用高斯窗口。我想引用起源,但现在找不到。该过程的思想在代码中给出。 看起来像这样:
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import cumtrapz
from scipy.optimize import curve_fit
import sys
sys.path.insert(0, "/home/nikousr/devel/2021-11-11_linearity_tester_python3/trunk/src/common/" )
from WaveAnalysis import Wave
def signal( x, A, B, C, D, W ):
r = (
A
+ B * x
+ C * np.sin( W * x )
+ D * np.cos( W * x )
)
return r
def signal_p( x, A, B, C, F, P ):
r = (
A
+ B * x
+ C * np.sin( F * x - P )
)
return r
testparams = [ -26, 2.03, 2.5, -1.56, 11.67]
### test data with noise
xl = np.linspace( -0.3, 2.6, 190)
sl = signal( xl, *testparams )
sl += np.random.normal( size=len( xl ), scale= 0.05 )
### numerical integrals
Sl = cumtrapz( sl, x=xl, initial=0 )
SSl = cumtrapz( Sl, x=xl, initial=0 )
### fitting the integro-differential equation to get the frequency
"""
note:
with y = A x**2 +...+ D sin() + E cos()
the double integral int( int (y ) ) = a x**4 + ... - y/F**2
"""
VMXT = np.array( [ xl**3, xl**2, xl, np.ones( len( xl ) ), sl ] )
VMX = VMXT.transpose()
A = np.dot( VMXT, VMX )
SV = np.dot( VMXT, SSl )
AI = np.linalg.inv( A )
result = np.dot( AI , SV )
print ( "Fit: ",result )
F = np.sqrt(-1/result[-1])
print("F = ", F)
### Fitting the linear parameters with the frequency known
VMXT = np.array(
[
np.ones( len( xl ) ),
xl,
np.sin( F * xl), np.cos( F * xl)
]
)
VMX = VMXT.transpose()
M = np.dot( VMXT, VMX )
SV = np.dot( VMXT, sl )
MI = np.linalg.inv( M )
A, B, C, D = np.dot( MI , SV )
print( " simple linear:")
print (A, B, C, D, F)
### Non-linear fit with initial guesses
### Try to make nonlinear fits with phases instead!
amp = np.sqrt( C**2 + D**2 )
phi = -np.arctan( C / D )
opt, cov = curve_fit( signal_p, xl, sl, p0=(A, B, amp, F, phi ) )
print( "non-linear fit:")
print( opt )
print("amplitudes of non-linear fit:")
print([opt[2] * np.cos( opt[-1]),opt[2] * np.sin( opt[-1])])
# ~ exit(0)
#### no slope no offset
sl2= sl - A - B * xl
dt = xl[1] - xl[0]
SR = 1 / dt
mywave = Wave( sl2, SR )
fr = mywave.get_base_frequency()
print(fr, 2 * np.pi * fr)
pssF, pssS = mywave.power_spectrum_single_sided()
### plotting
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot(
xl, sl,
ls='', marker='+', label="data", markersize=5
)
ax.plot(
xl, signal( xl, A, B, C, D, F ),
ls="--", label="double linear fit"
)
ax.plot(
xl, signal_p( xl, *opt),
ls=":", label="non-linear"
)
ax.legend( loc=0 )
ax.grid()
fig2 = plt.figure()
ax2 = fig2.add_subplot( 1, 1, 1 )
ax2.plot(
xl, sl2,
ls='', marker='+', label="data", markersize=5
)
ax2.grid()
fig3 = plt.figure()
ax3 = fig3.add_subplot( 1, 1, 1 )
ax3.plot(
pssF, pssS,
ls='', marker='+', label="data", markersize=5
)
ax3.grid()
ax3.axvline( x=fr)
plt.show()
并提供
Fit: [ 3.38686502e-01 -1.30013264e+01 -8.02977875e+00 -1.40654130e+00 -7.28674835e-03]
F = 11.71475241627821
simple linear:
-26.0235630892575 2.049006435336821 2.412893203776346 -1.6726851936757097 11.71475241627821
non-linear fit:
[-26.000498 2.02844789 2.9422975 11.67322791 0.55934161]
amplitudes of non-linear fit:
[2.4939050514194734, 1.5612662031036137]
f, w:
1.8572716933690985 11.66958221521727
正如人们所看到的,最后一种方法为频率提供了最佳结果。
Wave方法如下:
# coding=utf-8
"""
NMI 2018-06-27
class to analyse wave data
"""
import numpy as np
from scipy.signal import flattop, hann, gaussian ## fft window functions
from random import random
class Wave( object ):
"""
class that simplifies some fft procedures on periodic data
"""
def __init__( self, waveData, sampleRate ):
self._sampleRate = sampleRate
self._sampleSpacing = 1 / self._sampleRate
self._wave = np.array ( waveData, dtype=np.double ) ## real signal
self._wavePoints = len( waveData )
self._timelist = None
self._flist = None
self._fft = None
def wave( self ):
return( self._wave )
def number_of_points( self ):
return self._wavePoints
def sample_rate( self ):
return self._sampleRate
def sample_spacing( self ):
return self._sampleSpacing
def make_time_x( self ):
"""
generate an x-axis according to the given sample rate
"""
if self._timelist:
out = self._timelist
else:
data = np.arange( self._wavePoints )
out = data * self._sampleSpacing
self._timelist = out
return out
def make_freq_x( self ):
"""
provides the list of frequencies that correspnd to the fft
values follow the order (0, ...(n - 1 )/2, -( u ), ..., -1 ),
where u depends on whether n is even or odd.
see numpys fftfreq()
"""
if self._flist is None:
out = np.fft.fftfreq(
self._wavePoints,
self._sampleSpacing
)
self._flist = out
else:
out = self._flist
return out
def make_window( self, windowType=['FlatTop'], **kwargs):
"""
window functions for fft
"""
if windowType[0] == 'FlatTop':
window_func = flattop
elif windowType[0] == 'Hann':
window_func = hann
elif windowType[0] == 'Gaussian':
window_func = gaussian
else:
window_func = np.ones
window = window_func( self._wavePoints, *(windowType[1:]) )
normalisation = self._wavePoints / window.sum()
out = np.array( window ) * normalisation
return out
def fft( self, **kwargs ):
"""
fft of the wave using a window function given in kwargs.
Additional arguments are passed to the window function itself.
"""
out = np.fft.fft( self._wave * self.make_window( **kwargs ) )
return out / self._wavePoints
def power_spectrum_single_sided( self, **kwargs ):
powerS = self.fft( **kwargs )
powerS = ( powerS * np.conj( powerS ) ).real
powerS = 4 * powerS[ : self._wavePoints // 2 ]
powerF = self.make_freq_x()
powerF = powerF[ : self._wavePoints // 2 ]
return powerF, powerS
def spectrum_single_sided( self, **kwargs ):
pf, ps = self.power_spectrum_single_sided( **kwargs )
return pf, np.sqrt( ps )
def get_base_frequency( self, relativeSigmaQuotient=8 ):
"""
uses the fact that the fourier of a gaussian is a gaussian.
folding the virtual deltapeak of a frequency results in a gaussian.
taking the log and looking at two consecutive bins, preferably near the maximum,
allows to calculate the position of the maximum, i.e. the frequency
(square terms of the gaussians cancel and a linear equation is solved)
8 is my choice ... seems to work well...
just narrow enough to be zero on the edges and not mix peaks
"""
sigma = self._wavePoints / relativeSigmaQuotient
tau = self._wavePoints / ( 2.0 * np.pi * sigma )
pf, ps = self.spectrum_single_sided(
windowType=[ 'Gaussian', sigma ]
)
n = np.argmax( ps ) - 1
slg = tau**2 * np.log( ps[ n + 1 ] / ps[ n ] )
theoreticalbin = slg + n + 0.5
return theoreticalbin * self._sampleRate / self._wavePoints
评论