Estimation of Amplitude, Frequency and Phase of Linear Combination of Harmonic Signal Beyond the Leakage Resolution of DFT

Question

How can I find a rough ( as accurate as possible) Amplitude of each frequency when there is spectral leakage. Currently, I am dealing with a system that contains special leakage which seems unavoidable as I am measuring a real signal with any possible frequency. currently, I can identify the component frequencies to a decent level of accuracy even with the leakage, but I can't find the amplitude of this frequency as the energy is spread over multiple frequency bins.

Is it possible to get a rough prediction of a frequency amplitude either directly from the FFT, or afterwords using the identified frequencies?

I have been experimenting with Windowing and found this useful. https://www.youtube.com/watch?v=VxTx9QW8Zx8&ab_channel=Adash using the Hann window and using a frequency correction equation. it gives good approximations for frequency but no amplitude or phase.

A sample of the code currently used to test and generate the spectra.

import numpy as np

def gendata(FreqR,AmplR ,Ssize):
    PhaseR = np.pi # phase range
    FreqR = (FreqR*10) +1
    AmplR = (AmplR*100)+1
    PhaseR = (PhaseR*100)+1
    
    freq_OR = []
    ampl_OR= []
    phase_OR =[]


    for i in range(Ssize):
        freq_OR. append((((FreqR-1)/Ssize)*(i+1)/10)) # currently just equaly spread but can be any value between 0-9 rad/s
        ampl_OR.append(1) # unit amplitude to check for variation and error.
        phase_OR.append(np.random.randint((-PhaseR),(PhaseR))/100) # random phase 

    return(freq_OR,ampl_OR,phase_OR)


def pltsignal(freq_OR,ampl_OR,phase_OR,t_OR):
    
    z = 0 # original signal flat sea's
    if len(freq_OR) == len(ampl_OR):
        for i in range(len(freq_OR)): # for each frequency generate a regular wave.
            wave_OR = ampl_OR[i] * np.cos(2*np.pi*freq_OR[i]*t_OR + phase_OR[i]) # generated wave signal using the frequncy, amplitude and phase of each wave.
           
            z = z + wave_OR # superpostion each wave ontop of each other.
    else:
        print("amplutide and frequency are diffrent lengths")
 
    return z  # returns superimposed irregular wave.

FreqR = 1.5 # htz or (9 rad/s ish)
AmplR = 1 # currently unit amplitude to measure the error
Size = 12 # sample size, also used to move around frequencies

E_OR = 20 # time over wich signal is measured (max values is about 120 ish)
Fs = 6 # Sample rate in Hz

t_OR = np.arange(0,E_OR,1/Fs) # starts at 0, ends at E, steps by 1/Fs

data_OR = gendata(FreqR,AmplR,Size)
z = pltsignal(data_OR[0],data_OR[1],data_OR[2],t_OR) # original signal ( simulation of mesurment )

# anaysisng the signal using NumPy fft, then scaling and finding peaks.
Ramp = np.fft.fft(z) #real aplitudes used to find phase
Rfeq = np.fft.fftfreq(z.shape[-1]) #real frequency domain to find phase

Signal Model:

$$ x \left( t \right) = \sum_{i = 1}^{M} {a}_{i} \cos \left( 2 \pi {f}_{i} t + {\phi}_{i} \right) + n \left( t \right) $$

Where $ M, {\left\{ {a}_{i} \right\}}_{i = 1}^{M}, {\left\{ {f}_{i} \right\}}_{i = 1}^{M}, {\left\{ {\phi}_{i} \right\}}_{i = 1}^{M} $ are unknown parameters and $ n \left( t \right) $ is Additive Gaussian White Noise (AWGN).

One could assume:

1. The SNR is very very high.
2. The observation time is ~ 120 [Sec].
3. The number of signals is 1-20.
4. The frequencies are up to 2 [Hz].
5. The gap between frequencies can be as small as 0.005 [Hz].

Answer 1

Solving the Linear Combination of Real Harmonic Signal

The data model is given by:

$$ x \left( t \right) = \sum_{i = 1}^{M} {a}_{i} \sin \left( 2 \pi {f}_{i} t + {\phi}_{i} \right) + n \left( t \right) $$

Where $ M, {\left\{ {a}_{i} \right\}}_{i = 1}^{M}, {\left\{ {f}_{i} \right\}}_{i = 1}^{M}, {\left\{ {\phi}_{i} \right\}}_{i = 1}^{M} $ are unknown parameters and $ n \left( t \right) $ is Additive Gaussian White Noise (AWGN).

Estimating the parameters above requires 4 steps solution:

1. Estimate the model order ($ M $) / number of signals.
2. Estimate the frequencies accurately using direct method (Given the number of frequencies).
3. Estimate the amplitude and phase of the model given the frequencies (Usually using Non Linear Least Squares solver).
4. Estimate all amplitudes, frequencies and phase (Usually using Non Linear Least Squares solver).

The above works only in very very high SNR cases. For low SNR this problem is probably unsolvable (Certainly for large $ M $).

If the frequencies are well spread (Far from each other) then steps (1) and (2) can be achieved using DFT (fft()).
Yet in case the frequencies are close in the best case the estimation of the frequency won't be accurate (Leakage from near by frequencies will shift the peak) and in the worse the case frequencies will be blended and neither the number of signals nor the frequencies is achievable.

In the images above I used observation time of 12 [Sec] with 4 signals with frequencies of [0.95; 1.00; 1.05; 1.10] [Hz]. Since I took 1 / 10 of observation time in your question, it is equivalent of frequencies being 0.005 far apart in the 120 [Sec] observation time. I just wanted to deal with less samples.

As can be seen in the image above, in the DFT there is no way estimating the number of signals and certainly extracting the frequencies.

But using Super Resolution methods one could solve this case!
See some information about Frequency Domain Super Resolution in the section below.

Yet for this case I used different method for Super Resolution which I can't reveal (Used commercially by me as an advisor).
I also developed a method to estimate the number of signals.

In the case above indeed my method estimated 4 signals and the estimated frequencies were [0.9503, 1.0016, 1.0483, 1.0996] [Hz].

Then applying the next 2 steps I could even farther improve results:

As can be seen from above, the estimated signal is almost perfectly aligned with the ground truth (The model with known parameters with no noise).

Indeed the SNR here is very high yet still, without the Super Resolution method, it wouldn't work.

The full code is available on my StackExchange Signal Processing Q74024 GitHub Repository (Look at the SignalProcessing\Q74024 folder). Unfortunately, I can't disclose the code of the estimation of the number of parameters and the Super Resolution.

Frequency Domain Super Resolution

Old non modern methods will use some Windowing tricks which falls short in most cases. What you really need is Super Resolution. Super Resolution basically is equivalent to a "Window Method" with very narrow main lobe and almost no side lobes.

You can employ Compressed Sensing / Sparse Representation for Super Resolution in Frequency Domain.

One way to do so is solving the problem:

$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| F \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$

Super Resolution means, in that context, being able to resolve frequencies which are closer than what the observation time suggests (overcoming Leakage issues):

In the above you can see the DFT of a sum of 2 sines with the given frequencies. The Gaussian model is using $ {L}_{2} $ for regularization (Which is basically damped zero padding).

You may see that the $ {L}_{1} $ could resolve the 2 sines even when they are only 0.5 [Hz] apart with an observation windows of 1 [Sec].

The answer was taken from my answer for Super Resolution in Frequency Domain Using Compressed Sensing.

Answer 2

Applying ssq_cwt with extract_ridges, I obtain below. Improvable with better windowing, more samples. Smoothing can be applied on amplitude plot to make it more interpretable without losing much accuracy.

import numpy as np
from ssqueezepy import ssq_cwt, extract_ridges
from ssqueezepy.visuals import plot, imshow

# z = see OP's code; used np.random.seed(1000)

beta = 24
Tx, Wx, ssq_freqs, scales, *_ = ssq_cwt(z, ('gmw', {'beta': beta}), padtype='zero', fs=6)
ridge_idxs = extract_ridges(Tx, scales, penalty=20)

plot(ridge_idxs, color='k', linestyle='--', xlims=(0, len(z) - 1))
imshow(Tx, abs=1, yticks=ssq_freqs[::-1], ylabel="Frequencies [Hz]",
       title="abs(SSQ_CWT), wavelet=('gmw', {'beta': %s})" % beta)
amplitude = np.abs(Tx[ridge_idxs[:, 0], np.arange(len(z))])
frequencies = ssq_freqs[::-1][ridge_idxs[:, 0]]

plot(amplitude, ylims=(0, None), title="Amplitude vs time, SSQ ridge", show=1)
plot(frequencies, ylims=(0, None), ylabel="Frequencies [Hz]", 
     title="Frequency vs time, SSQ ridge", show=1)

Answer 3

FFT is poorly-equipped for this task$^{1}$; time-frequency localization, like STFT or CWT, is preferred. Said representations can be refined further to trace out frequency and amplitude over time with synchrosqueezing.

1: I originally understood the question as tracking instantaneous amplitude and frequency modulations of a (non-stationary) signal, which isn't the case; see below.

ssqueezepy offers all of the above, with wide variety of wavelet and windowing choices - ridge extraction included. Also see comparison of transforms and wavelet choices here - example:

---

UPDATE: OP's model is given by

$$ \sum_{i=0}^{M} a_i \cos(2\pi f_i t + \phi_i) $$

where $\phi_i$ is random. I.e., sum of fixed-amplitude, fixed-frequency sinusoids, with phase perturbations. Indeed FFT is finely equipped for this, and an "intelligent FFT" that circumvents resolution limitations will likely work better than time-frequency localization, as latter uses bandpass filters that suffer from classical resolution limitations.

(I can't speak with certainty as I don't understand "super resolution"; synchrosqueezing also 'improves' resolution its own way, and the most obvious advantage is redundancy that enables lot more powerful analysis methods.)

Answer 4

The standard method of dealing with spectral leakage is time domain windowing. This involves a fair bit of tradeoff: main lobe width, side lobe peaks & distribution, stop band attenuation, etc. These tradeoffs are controlled by choosing the window type and window parameters (if applicable).

What the best trade off is, really depends on your specific application requirements and types of signals.

Answer 5

Actually, there is a paper that attempts to addresses your question: Super-Resolving a Frequency Band The techniques involved may help you.