# Reconstructing a sine wave from an interval shorter than half its wavelength

I know that the function that I search is a sine wave of the form $$A \cdot \sin(x)$$

where my $$A$$ is unknown. I have samples of the sine wave at discrete points within an interval that's shorter than $$\lambda/2$$ represented by the red line in the plot. It must be possible to somehow determine $$A=1$$ through reconstruction or so from my samples given by the red line. How can i do that?

• Is the frequency of the sine wave known or is it also variable ? – Ben Aug 14 '20 at 14:00
• The frequency of the sine wave is known. – Bulbasaur Aug 14 '20 at 14:01
• Matlab code for the plot: f=1; w=2*pi*f; fs=1000; A=1; t=0:1/fs:1; func=A*sin(w*t); plot(func) title('sine wave') ylabel('amplitude') xlabel('sample') hold on t2=(0:1/fs:0.2); func2=A*sin(w*t2); plot(func2) – Bulbasaur Aug 14 '20 at 14:13
• en.wikipedia.org/wiki/Sinusoidal_model – endolith Aug 14 '20 at 14:13
• My intuition is that you could try a least-squares fit on that curve. However, I don't think it's gonna be robust to noise, harmonics, DC offset, etc. – Ben Aug 14 '20 at 14:22

import numpy as np

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

omega = 0.1

X = np.array( [ 0.5, 0.6, 0.65, 0.7 ] )

C = np.zeros( 4 );  S = np.zeros( 4 );  U = np.zeros( 4 )

for n in range( 4 ):
C[n] = np.cos( omega * n )
S[n] = np.sin( omega * n )
U[n] = 1.0

V = np.zeros( 3 )
M = np.zeros( ( 3, 3 ) )

M[0,0] = C.dot( C );  M[0,1] = S.dot( C ); M[0,2] = U.dot( C )
M[1,0] = M[0,1];      M[1,1] = S.dot( S ); M[1,2] = U.dot( S )
M[2,0] = M[0,2];      M[2,1] = M[1,2];     M[2,2] = U.dot( U )

V = X.dot( C );    V = X.dot( S );   V = X.dot( U )

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

for n in range( 4 ):
y = R * C[n] + R * S[n] + R * U[n]
print ( n, X[n], y )

#========================================================================
main()
0 0.5 0.502540182211
1 0.6 0.592404834029
2 0.65 0.657595165971
3 0.7 0.697459817789

• I have a follow up question to this script. I implemented it in Matlab and it works great for low frequencies. But as soon as I increase the frequency above 100 Hz, the entries in C become about 1 and those in S become about 0. The matrix M becomes singular and the entries in R are infinite. I dont quite understand whats happening. Do you have an idea why this problem occurs and how I can solve that problem? – Bulbasaur Aug 31 '20 at 7:58
• If you are close to a peak, and take a very short stretch, then your cosine basis (close to all 1s) is not going to be independent of your DC (all 1s), and hence the problem. Simply drop the DC basis vector. This forces your solution to assume it is zero based and you lose the DC offset value. – Cedron Dawg Aug 31 '20 at 12:23

If your signal is really as simple as

$$x(t)=A\sin(\omega_0t)\tag{1}$$

with known $$\omega_0$$, and you have observations $$y(t_i)$$, which are noisy samples of $$x(t)$$ at known time instances $$t_i$$, then a simple solution would be the least squares estimate

$$\hat{A}=\frac{\displaystyle\sum_iy(t_i)\sin(\omega_0t_i)}{\displaystyle\sum_i\sin^2(\omega_0t_i)}\tag{2}$$

Of course, this simple solution won't work if your signal actually has the form

$$x(t)=A\sin(\omega_0t+\phi)+c \tag{3}$$

with unknown phase $$\phi$$ and DC-offset $$c$$. However, you can compute optimum least squares estimates for that problem too. This is discussed in Cedron's answer.

• Actually the samples are not even noisy in my case. The least square estimate works suprisingly well (I have no experience with it so far but it returns exactly A=1. Thanks. – Bulbasaur Aug 14 '20 at 16:23
• @Bulbasaur: If you have no noise and no other errors of any kind, you don't even need the sums, you might as well just take any value $t_i$ and compute $A=y(t_i)/\sin(\omega_0t_i)$. Of course, this only works as long as you know the frequency $\omega_0$ and the phase exactly. However, I don't think that this is a realistic problem. – Matt L. Aug 14 '20 at 16:26
• In my special case the only "Error" is machine precision, so I will probably still be a little better of with the sum. Your answer certainly has been helpful to me. – Bulbasaur Aug 14 '20 at 16:30
• @Bulbasaur Matt's simpler answer gives your point at (0,0) superpowers and can't be budged. Mine calibrates for any offsets and will tell you where it thinks (0,0) is. The relative sizes of the coefficient values should reflect this, a,c << b. If your samples were only near the peak, this could make the two alternatives noise sensitive and if your set says (0,) is here, so to speak, you would want that anchor. – Cedron Dawg Aug 14 '20 at 17:01
• @Bulbasaur: Yes, that's clearly the case. But you should add the exact specification to your question instead of disclosing it bit by bit :) – Matt L. Aug 14 '20 at 17:37

Build a basis set with your frequency and match your signal. It is straightforward linear algebra:

$$C$$ is portion of cosine

$$S$$ is portion of the sine

$$U$$ is a vector of ones (DC)

$$X = a C + b S + c U$$

$$X \cdot C = a (C \cdot C) + b (S \cdot C) + c (U \cdot C)$$ $$X \cdot S = a (C \cdot S) + b (S \cdot S) + c (U \cdot S)$$ $$X \cdot U = a (C \cdot U) + b (S \cdot U) + c (U \cdot U)$$

Now you have three equations three unknowns, $$a$$, $$b$$, and $$c$$.

Best fit interpolation/extrapolation function:

$$x[n] = a \cos[wn] + b \sin[wn] + c$$

$$A = \sqrt{a^2+b^2}$$

Now, wouldn't it be handy if $$C\cdot S=0$$?

[Overengineered solution to account for any vertical or horizontal shifts, use Matt's if you know it is a simple multiple.]

• won't you get an overdetermined set of equations? – Ben Aug 14 '20 at 15:38
• @Ben You got your starting X vector, you create a C and S (which happen to match the frequency of your target) and you are doing a standard least squares fit. (Matt's answer is $b = X\cdot U / S\cdot U$ where $U$ is all ones, aka DC of one.) U would be the first vector to use and should be included. Oops, I'll correct that. – Cedron Dawg Aug 14 '20 at 16:04
• @Ben: The first equation $X=aC+bS+cU$ is the overdetermined system you're looking for. – Matt L. Aug 14 '20 at 19:15
• @MattL. Overdetermined? Everybody should see this as a proto-DFT. – Cedron Dawg Aug 14 '20 at 20:09
• @Ben Matt changed his answer. It is now $b = X\cdot S / S \cdot S$. Shouldn't change the answer any, but extra credit, what is the noise impact? – Cedron Dawg Aug 14 '20 at 22:59

If samples are results of precise measurements, the amplitude is any of ratios $$sample_i/sin(ω·t_i)$$.

If the noise is present, some kind of averaging is needed. Because a single parameter is calculated, and no information is given on noise distribution, a simple weighted average is the only available option: $$A_{avg} = {Σ(sample_i·sin(ωt_i))\over{Σsin^2(ωt_i)}}$$. A weighted sample variance is $$(σ_w)^2 = {Σ\{sin(ωt_i)·(sample_i-A_{avg}·sin(ωt_i))^2\}\over{Σsin(ωt_i)}}$$

• you need to know the amplitude of the sine wave to calculate the arcsine of each sample... – Ben Aug 14 '20 at 14:52
• Yes, it's been a flop of a kind. Hope this correction will do.. – V.V.T Aug 14 '20 at 16:17