[time-nuts] simple phase finder

[Detlef.Amberg at gmx.de]

Hi,

if you need the frequency, amplitude and phase of a sampled sine wave 
there is a solution without correlation, nonlinear fit or interpolation 
in the frequency domain. It works with damped sine waves as well.

You just have to solve a linear system of equations.

See https://www.dsprelated.com/showarticle/795.php
with links to Matlab code.

Cheers
Detlef


Am 05.12.2018 um 14:19 schrieb jimlux:
> 
> Here's some python code that I use to find the frequency, amplitude and 
> phase of a discrete sine in noise for 100MHz samples. It's basically a 
> sequential sliding correlator.
> 
> If I were decoding WWVB to start, I'd break my samples up into 0.1 
> second or 0.5 second chunks and process them to see what the carrier 
> phase is.  If I did 0.1 second chunks, I can probably identify the bit 
> transitions every second, because in 1/10th chunks, the phase will not 
> be one of two values.
> 
> This routine is not computationally efficient, it's pretty brute force - 
> a better approach would use a narrow band transform and do a correlation.
> 
> This routine also might get fouled up by the amplitude modulation (at 
> least in the "peak search" part of operation.
> 
> Another approach would be to implement a classical Costas Loop.  I think 
> though, a sliding correlator of some sort might be a better solution - 
> for one thing, it "look forward and back in time"
> 
> 
> fs = sample rate
> M = number of samples
> ftestmin, ftestmax is the range to search over in MHz
> adc is a numpy array with the samples
> 
> 
> 
>      # build an array of sample times
>      t = np.arange(0,M)
>      t = t/fs
>      # iteratively search a small range around the peak to find the best 
> fit for a sine wave.
>      # the resolution bandwidth for 32768 points is about 3 kHz, so 
> looking over
>      # to make life nicer, we'll round the start and stop frequency to a
>      # multiple of 100 Hz, then go in 10 Hz steps
>      ftestmin = 0.0001 * math.floor(ftestmin * 10000)
>      ftestmax = 0.0001 * math.ceil(ftestmax * 10000)
> 
>      resid = adc - np.mean(adc)
> 
>      ftest = np.arange(ftestmin,ftestmax,0.000010)
>      test1   = 0
>      testmax = 0
>      pi = np.pi
>      for i in range(0,len(ftest)) :
>          try1 = (np.cos(t * 2 * pi * ftest[i]) - 1.0j*np.sin(t * 2 * pi 
> * ftest[i]))
>          try1 = np.reshape(try1, (adc.size,1))
>          test1 = np.sum(resid * try1,axis=0) / M
>          if abs(test1[0]) > abs(testmax):
>              testmax = test1
>              ftestmax = ftest[i]
> #
> 
>      c = np.cos(t * 2 * pi * ftestmax)
>      s = np.sin(t * 2 * pi * ftestmax)
> 
>      f1db = 20 * np.log10(np.sqrt(2) * abs(testmax))
> 
> 
>      freqreturn  = ftestmax
>      ampreturn   = f1db[0]
>      phasereturn = np.angle(testmax[0]) * 180 / pi
> 
> _______________________________________________
> time-nuts mailing list -- time-nuts at lists.febo.com
> To unsubscribe, go to 
> http://lists.febo.com/mailman/listinfo/time-nuts_lists.febo.com
> and follow the instructions there.
>