我们如何准确估计这种信号的频率?

How can we accurately estimate the frequency of such a signal?

提问人:BlueConda 提问时间:8/7/2023 更新时间:8/31/2023 访问量:75

问:

我有一个信号,其表达式如下所示:

信号的数学表达式

当信号中的ABωDE都是未知参数时,我可以使用什么方法来准确拟合该信号?特别是,其中最关键的未知参数是ω。

信号如下所示:

信号是什么样子的

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

我尝试了 MATLAB 的 lsqcurvefit 和 nlinfit 函数,但准确性很难令人满意。尝试使用 FFT 提取 ω,发现结果也不是很准确。

Python MATLAB 曲线拟合 频率分析

评论

1赞 mozway 8/7/2023
您的真实数据是怎样的?请提供一个可重复的例子。拟合的准确性很大程度上取决于数据的噪声和采样率。
0赞 Irreducible 8/7/2023
嗨,欢迎来到 SO,请阅读 如何提问最小重现示例。只是说你所有的方法都不令人满意,而不展示你做了什么来估计参数,这对我们没有帮助,对你也没有帮助。请提供所有必要的信息。
1赞 Luis Mendo 8/8/2023
也许高通滤波器去除了线性和常数项,然后你可以更可靠地使用 FFT 来提取 omega?要应用高通滤波器,您需要知道 omega 的一些下限
0赞 BlueConda 10/29/2023
非常感谢您的热情回复。你的回答给了我很好的解决方案。我非常感激。

答:

1赞 John BG 8/27/2023 #1

1.- 问题所附图中提供的信号清楚地显示 单音 + 1 锯齿。

这是提供的“屏幕截图”

enter image description here

只要音调有足够的功率显示在锯齿波频谱以上 然后估计单个音调的寻求频率,只需发现 被锯齿光谱包围的单音的峰值。

如果所寻求的音调频率发生在功率低于锯齿波的包络线时,则必须沿着锯齿波谱包络线寻找异常,否则,如果没有音调,则随着频率的增加而平滑衰减。1/k

2.- 您正在分析的信号与提供的信号不完全相同 表达

这是提供的信号表达式

enter image description here

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

改写;信号没有傅里叶变换,或者因为变换到频谱的积分不产生任何东西。x1(t)=a*tx2(t)=a*|t|[-Inf Inf][0 Inf]

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

由于没有提到或时间参考,我假设一个:tstarttstopt

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*tstopmin(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 个锯齿。

enter image description here

enter image description here

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]');

enter image description here

enter image description here

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]');

enter image description here

enter image description here

因此,很明显,问题中提供的信号图像是 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');

enter image description here

enter image description here

与单音相比,锯齿波具有填充的频谱,除了噪声和干扰之外

enter image description here

enter image description here

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')

enter image description here

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

例如,我们可以回顾锯齿信号的基础知识:

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

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

6.- 锯齿信号频谱包络

T2是锯齿基循环。

如果锯齿被时间定义为无直流电,则删除直流常数。

A2是我最初定义的带直流的锯齿振幅。

零直流锯齿,具有振幅。A2/2

锯齿信号(非功率)频谱的包络是,但这不包括和:1/(pi*k)FsL

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)

enter image description here

enter image description here

评论

请注意,没有 DC,因为音调是,我已经移除了任何锯齿 DC。sin

对于 2 个音调,由于两者都在同一频率上,因此幅度谱相似。

如果您想了解如何捕捉带有锯齿信号的微弱音调,请发布另一个问题。

0赞 mikuszefski 8/31/2023 #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

Noisy data with fits FFT with best frequency

正如人们所看到的,最后一种方法为频率提供了最佳结果。

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