[time-nuts] simple phase finder
Detlef.Amberg at gmx.de
Detlef.Amberg at gmx.de
Wed Dec 5 23:11:30 UTC 2018
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.
>
More information about the Time-nuts_lists.febo.com
mailing list