[time-nuts] simple phase finder

[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