2
$\begingroup$

Before you read my problem, I've summed up my experiments in a notebook available here.

So, I have a few noisy samples (~50) of a pseudo sinusoidal signal for which I use a traditional discret Fourier transform in order to extract the amplitude and phase of the fundamental (real signal is much noisier).

fft[1]=5.0, 60.0°

enter image description here Now my problem is that i sometime have only half the data available, ie only half a peridod of that sinusoid. But I know I almost exactly have sample of half that period, and I know the signal is almost a sinusoid.

How can I proceed to extract the fundamental characteristics anyway?

I've noticed that zero-padding the missing half and doubling the amplitude works quite well.

fft[1]=2.5, 59.1°

enter image description here

But in my case the sinusoid is not DC-centered on zero, and in that case zero-padding doesn't work.

fft[1]=10.6, -82.5°

enter image description here

I would need to pad with the DC of my signal, which is unknown. I've noticed that if you have half a period, the average of first and last samples should give the DC, but this won't be reliable with noisy data.

fft[1]=2.5, 59.1°

enter image description here

Thanks in advance for your insights!

$\endgroup$
3
$\begingroup$

If you know the frequency (and it seems that you do), simply construct a basis set of vectors from

1) DC (all ones)

2) Half a Sine wave

3) Half a Cosine wave

These form an orthogonal basis just like a DFT. In fact you can think of it as a "one half bin and DC bin" specialized DFT.

Discard the DC reading (unless you care) and use the Sine and Cosine coefficients just like a regular DFT bin from which you can calculate magnitude (which gives you amplitude with a multiplication) and the phase (using atan2).

This will give you a best fit answer. No window function and no zero padding.


Small correction (sorry): The basis is not orthogonal. You will need to solve a system of equations.

$$ x[n] = A \cos( \omega n ) + B \sin( \omega n ) + C $$

As a vector equation:

$$ \vec x = A \vec c + B \vec s + C \vec u $$

With some dot products, you get the system of equations with $A,B,C$ as the unknowns to be solved for.

$$ \vec x \cdot \vec c = A \vec c\cdot \vec c + B \vec s\cdot \vec c + C \vec u \cdot \vec c $$

$$ \vec x \cdot \vec s = A \vec c\cdot \vec s + B \vec s\cdot \vec s + C \vec u \cdot \vec s $$

$$ \vec x \cdot \vec u = A \vec c\cdot \vec u + B \vec s\cdot \vec u + C \vec u \cdot \vec u $$

Note:

$$ \vec c \cdot \vec u = \vec u \cdot \vec c = 0 $$

$$ \vec c \cdot \vec s = \vec s \cdot \vec c = 0 $$

$$ \vec u \cdot \vec u = N $$

Where $N$ is your sample count. Thus the equations simplify to:

$$ \vec x \cdot \vec c = A \vec c\cdot \vec c $$

$$ \vec x \cdot \vec s = B \vec s\cdot \vec s + C \vec u \cdot \vec s $$

$$ \vec x \cdot \vec u = B \vec s\cdot \vec u + C N $$


In order to get the best results with noisy data it is best to build the basis vectors over the full range of data you have and then solve the full set of equations. A little more work, but likely worth it.


Okay, for the OP and maybe Ben. Here is a sample program for the best fit of a partial cycle. In this case 0.6 of a cycle.

import numpy as np
import matplotlib.pyplot as plt

#==========================================================
def main():

#---- Set Parameters, Base Frequency and Phase

        f_s = 50.0
        p   = 60.0 * np.pi / 180.0

#---- Set Domain Time and Angle for partial cycle

        t = np.arange( 0.0, 1.0, 1.0 / f_s )
        a = 2.0 * np.pi * t * 0.6

#---- Construct the Signal

        x = 10.0 * np.cos(        a + p ) \
          +  0.2 * np.sin(  2.0 * a )     \
          +  0.2 * np.sin( 10.0 * a )

#---- Build Basis Vectors

        c = np.cos( a )
        s = np.sin( a )
        u = np.ones( a.shape )

#---- Do Basis Dot Products

        cc = c.dot( c )
        ss = s.dot( s )
        uu = u.dot( u )
        cs = c.dot( s )
        cu = c.dot( u )
        su = s.dot( u )

        print cs, ss, uu, cs, cu, su

#---- Do Signal Dot Products

        xc = x.dot( c )
        xs = x.dot( s )
        xu = x.dot( u )

        print xc, xs, xu

#---- Build Coefficient Matrix

        M = np.array( [[ cc, cs, cu ],  \
                       [ cs, ss, su ],  \
                       [ cu, su, uu ] ] )

        print M                       

#---- Build Values Matrix

        V = np.array( [xc, xs, xu] )

        print V

#---- Solve for the best fit "Bin Value"                       

        R = np.linalg.solve( M, V )

        print R

#---- Calculate the Best Fit parameters

        A = R[0]
        B = R[1]
        C = R[2]

        D   = np.sqrt( A*A + B*B )
        phi = np.arctan2( -B, A )

        print D, phi, p

#---- Construct the Best Fit signal

        y = D * np.cos( a + phi ) + C

#---- Show the Results

        plt.plot( t, x )
        plt.plot( t, y )

        plt.show()


#==========================================================
main()

And these are the results:

enter image description here

These are the values:

Underlying target:

    10.0
p =  1.0471975512

Found best fit:

D   = 10.1939671316 
phi =  1.04122684844 

Changing the parameters a little bit for a smaller portion (0.3) of the cycle and moving it up by 20.

        x = 10.0 * np.cos(        a + p ) \
          +  0.2 * np.sin(  2.0 * a )     \
          +  2.0 * np.sin( 10.0 * a )     \
          + 20.0

Results:

enter image description here

Notice this isn't as accurate because of the smaller cycle portion.

D   = 10.4409622651
phi =  0.943824736628

However, if the signal is a pure tone:

        x = 10.0 * np.cos( a + p ) + 20.0

The results are still exact on less than a third of the cycle:

D   = 10.0 
phi =  1.0471975512
| improve this answer | |
$\endgroup$
  • $\begingroup$ I'm not the OP, but I will try it $\endgroup$ – Ben Apr 21 at 1:10
  • $\begingroup$ @Ben Very good. I kind of blew the simplicity of the answer. See my update. $\endgroup$ – Cedron Dawg Apr 21 at 2:00
  • 1
    $\begingroup$ Ced, are you sure that $\vec u \cdot \vec s = 0$? i agree with you that $ \vec s \cdot \vec c = 0$ and $ \vec u \cdot \vec c = 0$. but i don't think that $ \vec u \cdot \vec s = 0$ for the half period that the phase is between $0$ and $\pi$. $\endgroup$ – robert bristow-johnson Apr 21 at 2:36
  • $\begingroup$ @robertbristow-johnson Did you add this comment before my followup? What you say is correct. $\endgroup$ – Cedron Dawg Apr 21 at 10:47
  • $\begingroup$ How do you justify the decomposition in $ A \cos( \omega n ) + B \sin( \omega n ) + C $, by using $ D \cos( \omega n + p) = D \cos( \omega n) \cos(p) - D \sin( \omega n) \sin(p) $ ? In that case $ A = D \cos(p), B = -D\sin(p) $ $\endgroup$ – Julien M Apr 22 at 7:15
1
$\begingroup$

Let's see if this algorithm is helpful to you.

  1. You could denoise the half period of samples you have available and then calculate the DC on this sinusoid using this denoised signal.

  2. Next you zero pad and interpolate the original half period signal (noisy version) and now pad the DC you calculated is step 1

| improve this answer | |
$\endgroup$
  • $\begingroup$ Can you tell me how to compute the DC in step 1, considering I have half a sinusoid with unknown phase? $\endgroup$ – Julien M Apr 19 at 7:17
  • $\begingroup$ So first step is to low pass filter the noisy signal, you could use a moving average filter or leaky integrator, any low pass filter that is good enough for the nosie you have, then integrate over half the period the modulus of the sinuoid/your filtered signal (since you do not know the phase), I assume you know the frequency, this result should always be constant over any phase, and the right hand side of this integral will depend on the amplitude (which would give you the DC, you can now calculate RMS) and frqeuency $\endgroup$ – Dsp guy sam Apr 19 at 7:31
  • $\begingroup$ Sorry to be slow. Yes I do know the frequency. Do you mean that after filtering, I should sum abs(signal) over 0-Pi? It's not clear how I can retrieve the DC. Integral of a sinusoid over half a period is not constant. And I can't take the modulus of the signal if I don't know the DC $\endgroup$ – Julien M Apr 19 at 16:12
  • $\begingroup$ I mean after denoising take the absolute value of the signal and integrate this function over half a period $\endgroup$ – Dsp guy sam Apr 19 at 16:19
  • $\begingroup$ My signal is strictly positive (values are pixel coordinates), if I sum its values over half a period, I won't obtain something constant over any phase. For instance ⌒ or ◡ won't sum up to the same value $\endgroup$ – Julien M Apr 19 at 21:15
0
$\begingroup$

Assuming that:

  • The input is 1/2 period, frequency is known
  • The waveform is exactly a sinusoid
  • The signal is sufficiently over-sampled
  • The amplitude is unknown
  • The phase is unknown

I guess you could do time-domain cross-correlation, something like this:

fs = 30;
T = 1;
t = 0:(1/fs):(T-1/fs);
a = 2.4;
phi = +1.1;
f = 0.5;
x = a*sin(2*pi*f*t+phi);
x = x + 0.1*randn(size(x));

t_template = 0:(1/fs):(T/((2/3)*f)-1/fs);
x_template = sin(2*pi*f*t_template);
[c, lags] = xcorr(x, x_template);
[mv, mvi] = max(c);
mvi = lags(mvi);
a_hat = sum(x.^2) / mv
phi_hat = -(pi/2)*(mvi/fs)/f
| improve this answer | |
$\endgroup$
  • $\begingroup$ I believe it won't work with my waveform which is a waveform+DC, DC contribution to xcorr will bias the max $\endgroup$ – Julien M Apr 22 at 7:07
  • $\begingroup$ Are you adding DC independently, or are you talking about the DC offset that stems from truncating a sinoid? $\endgroup$ – Knut Inge Apr 22 at 14:54
  • $\begingroup$ Independant (unknown) DC $\endgroup$ – Julien M Apr 22 at 22:23
0
$\begingroup$

This sounds like an instantaneous frequency estimation question. Have you ever tried Hilbert Transform? Or Hilbert Huang transform (For noisy signal )?

In fact, according to statistic signal processing, the CRLB of single tone frequency estimation problems hast been discussed and bounded via observing time (dominant), SNR, and fps. During the development of MLE, DTFT hast been taken into account so the zero padding (sampling more from DTFT) has its limits.

So, some no classical method (maybe parameter based ) should be tried.

Here I listed some method I know:

  1. AR/ARMA model

  2. Hilbert transform

| improve this answer | |
$\endgroup$
  • $\begingroup$ I quickly tried pyhht.emd.EMD() but it fails with an error since, with only half a period, it cannot find any IMF. Regarding AR/ARMA, i'm not sure it will help since I'm more interested in the magnitude/phase of my sinusoid than it's frequency. (I assume the frequency is known) $\endgroup$ – Julien M Apr 20 at 23:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.