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 -- [email protected]
To unsubscribe, go to
http://lists.febo.com/mailman/listinfo/time-nuts_lists.febo.com
and follow the instructions there.
