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

Question

I have a signal whose expression is shown below:

Mathematical expression of the signal

What method can I use to fit this signal accurately when A, B, ω, D,E in the signal are all unknown parameters? In particular, the most critical unknown parameter in there is ω.

The signal looks something like this:

What the signal looks like

How can the value of ω, the most critical of these, be estimated accurately?

I tried MATLAB's lsqcurvefit and nlinfit functions, but the accuracy is hardly satisfactory. Tried extracting ω using FFT and found that the results are not very accurate either.

Answer 1

1.- The signal supplied in the figure attached to the question clearly shows a single tone + 1 sawtooth.

This is the supplied kind of a 'screen shot'

As long as the tone has enough power to show above the sawtooth spectrum then to estimate the sought frequency of a single tone sufices to just spot the peak of the single tone surrounded by the sawtooth spectrum.

If the sought tone frequency happens with power below the envelope of the sawtooth, then you have to look for the anomaly along the sawtooth spectrum envelope that otherwise, without a tone, has a smooth 1/k decay with increasing frequency.

2.- The signal you are analysing is not exactly the one in the supplied expression

this is the supplied expression of the signal

Even if the supplied signal is all you have, one must assume finite power signals, not infinite energy signals, in order to use the Fourier transform.

Rewording; the signals x1(t)=a*t, x2(t)=a*|t| do not have a Fourier transform either along [-Inf Inf] or [0 Inf] because the integral for the transform to spectrum does not produce anything.

So we need to time window the signal and therefore the signal to be considered for spectrum analysis is sum of 1 tone and 1 sawtooth.

Since no tstart tstop or t time reference mentioned I asume one :

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

If f2*tstop not an integer then either a fraction of a min(f1,f2) should be added or cut last cycle fraction.

For simplicity I have not included such case where there are cycle tails to deal with.

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

So it's clear that the image of the signal supplied in this question is 1 tone and 1 sawtooth.

3.- The supplied image of the signal IS NOT 2 tones and 1 sawtooth

A signal including 2 tones both on same frequency and 1 sawtooth actually looks as follows :

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

So it's clear that the image of the signal supplied in the question is 1 tone and 1 sawtooth.

The time discontinuities, as time signals have been generated, can be solved so there are smooth transitions only. It's not difficult, but laborious, and for the purpose of answering this question, such effort would not add relevant contents.

4.- Spectrum

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

In contrast to single tones, sawtooths have populated spectrums, noise and interference aside

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.- References to Sawtooth signal basics:

We can review sawtooth signal basics in for instance:

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

and

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

6.- Sawtooth signal spectrum envelope

T2 is the sawtooth base cycle.

Remove the DC constant if the sawtooth was been time defined as with no DC.

A2 is the sawtooth amplitude WITH DC as I have initially defined it.

The null DC sawtooth, has amplitude A2/2.

The envelope of the sawtooth SIGNAL (NOT POWER) spectrum is 1/(pi*k) but this does not include Fs and 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)

Comments

Note there's no DC because the tone is sin and I have removed any sawtooth DC.

For 2 tones, since both would be on the same frequency the amplitude spectrum would be similar.

If you would like to see how to catch weak tones cluttered with a sawtooth signal please post another question.

Answer 2

Following the ideas of SE member Jean Jacquelin it easy to get a first estimate, assuming the data is dense.

In this case a double numerical integral of you function y, i.e. SSy, has the property SSy = -w^2 y + D/6 t^3 + E/2 t^2 + U t + Z. So with simple linear equations we get w. With known w the entire problem becomes purely linear and we can fit for the other parameters. The result can be improved by plugging this into a non-linear fit. One has to be careful, though. Sine-Cosine fits work better when introducing a phase instead ( see some details here ). Finally, we can try to improve on the frequency using FT. First I remove the offset and slope. For precise determination of the frequency I use a Gaussian window. I'd like to cite the origin, but can't find it right now. The idea of the procedure is given in the code. Looks like this:

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

and provides

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

As one can see, the last method provides the best result for the frequency.

The Wave-methods are as follows:

# 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