jimlux

2018-12-05 13:19:35 UTC

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-***@lists.febo.com

To unsubscribe, go to http://lists.febo.com/mailman/listinfo/time-nuts_lists.febo.com

and follow the instructions there.

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-***@lists.febo.com

To unsubscribe, go to http://lists.febo.com/mailman/listinfo/time-nuts_lists.febo.com

and follow the instructions there.