# Sine Wave autocorrelation, interpolation of phase

Started by July 31, 2017
Hallo.
For the purpose of measuring complex impedances i need to compare the
phase of two copies of a sinewave over a number of periods.
The sinewave is generated in the same &micro;C. Frequency is known and stable.
The original and the shifted signals are sampled by a double synchr. ADC.
The phase shift is the base for calculation of the complex Z of a load.
I have implemented an (auto-)correlation of the two signals. Works fine so
far.
Since my samplig-F is not so much higher than the sinewave, I need to
implement interpolation to get a higher resolution.
The autocorrelation of a sine wave is a cosine. Steepness around 0&deg; is
very low. What algo is recommended for interpolation?
Shall i use a different approach in general?
Thanks.
Eduardo

---------------------------------------
Posted through http://www.DSPRelated.com

On Monday, July 31, 2017 at 7:07:01 AM UTC-7, eduardoG26 wrote:
> Hallo.
> For the purpose of measuring complex impedances i need to compare the
> phase of two copies of a sinewave over a number of periods.
> The sinewave is generated in the same &micro;C. Frequency is known and stable.
> The original and the shifted signals are sampled by a double synchr. ADC.
> The phase shift is the base for calculation of the complex Z of a load.
> I have implemented an (auto-)correlation of the two signals. Works fine so
> far.
> Since my samplig-F is not so much higher than the sinewave, I need to
> implement interpolation to get a higher resolution.
> The autocorrelation of a sine wave is a cosine. Steepness around 0&deg; is
> very low. What algo is recommended for interpolation?
> Shall i use a different approach in general?

If you have enough samples, it should work, but it does depend
a little on what the sine actually looks like.

Is the complex load perfectly linear, the output will be a perfect
sine with phase shift and amplitude change.  But many loads aren't
quite linear enough for that.

One possibility is to do a non-linear least squares fit to
A*sin(B*t+C)  If B is well enough known, and since you are
generating the sine it should be, you could do a linear least
squares fit to D*sin(w*t) + E*cos(w*t).  (That is, linear in
the coefficients D and E.)  You can do this even if the sampling
rate is very low, if you have enough samples.

You could also add F*sin(3*w*t) + G*cos(3*w*t) terms as a
measure of the non-linearity.  Well, it depends on what
you are measuring, but many loads generate only odd harmonics.


If you only need to compare the phase, you could use the Hilbert transform on each of the signals. Extract the instantaneous phase from each signal. Subtract the phase and you will get the phase difference between the 2.


On Monday, July 31, 2017 at 10:07:01 AM UTC-4, eduardoG26 wrote:
>
> For the purpose of measuring complex impedances i need to compare the
> phase of two copies of a sinewave over a number of periods.
> The sinewave is generated in the same &micro;C. Frequency is known and stable.
> The original and the shifted signals are sampled by a double synchr. ADC.
> The phase shift is the base for calculation of the complex Z of a load.
> I have implemented an (auto-)correlation of the two signals. Works fine so
> far.
> Since my samplig-F is not so much higher than the sinewave, I need to
> implement interpolation to get a higher resolution.
> The autocorrelation of a sine wave is a cosine. Steepness around 0&deg; is
> very low. What algo is recommended for interpolation?
> Shall i use a different approach in general?

your sample rate still needs to be faster than twice that of the sinusoid.

you don't need to interpolate cross-correlation.  you need to evaluate the cross-correlation at zero lag and the AC powers of both the input and output signals.  and if you know approximately the impedance (such as positive real part and either capacitive or inductive imaginary part), then you can calculate both the magnitude and the precise phase from those three numbers.

Z =  |Z| e^{j \theta}

where |Z| = sqrt( R_vv(0) / R_ii(0) )

cos(theta) = R_vi(0) / sqrt( R_vv(0) R_ii(0) )

where R_vi(0) is the cross-correlation (at lag 0) of the voltage signal against the current signal, R_vv(0) is the autocorrelation of the voltage signal and R_ii(0) is the autocorrelation of the current signal.

make sure there is no DC in any of these signals.  there could be something, just from an offset in the A/D converters. use a DC-blocking algorithm on both v[n] and i[n], to make sure there is no DC.

r b-j

Thanks for the first three replies. Very helpful.
I will take a look at Hilbert which is new to me. Maybe a ready-to-go
solution. Thanks for the tip.
Resuming all three replies I will try the following.
Correlate both sine wave with sine and cosine (See I/Q demod.). Calc exact
phase by tangens(I/Q). Subtract phases. Ready. (Will this work?)

Note:
Sampling is 4Ms. Sine is ~40KHz. Will have to remove DC by SW. A great
advantage is the exact knowledge of sine period length in samples.
Components are linear. Combinations are R, C, L, R+C and R+L series and
parallel. Extracting the equations = a challenge. Got only RC series so
far. Still working on the other 3. Looks easier as I thought so far...

---------------------------------------
Posted through http://www.DSPRelated.com

On 31.07.17 16.06, eduardoG26 wrote:
> For the purpose of measuring complex impedances i need to compare the
> phase of two copies of a sinewave over a number of periods.
> The sinewave is generated in the same &micro;C. Frequency is known and stable.
> The original and the shifted signals are sampled by a double synchr. ADC.
> The phase shift is the base for calculation of the complex Z of a load.

Your setup is perfect. You can get phase information several orders of
magnitude below the sampling rate. BTDT. (I measured the length of an
audio cable with an ordinary 48 kHz on board sound device. Sampling rate
~21&#2013266101;s, time shift resolution <5ns)

> I have implemented an (auto-)correlation of the two signals. Works fine so
> far.

But I would not use sinusoids for this purpose. Measuring at exactly one
frequency lets you get either only unspecific, noisy information at the
ground frequency of your FFT length or, at higher frequencies, many
ambiguous results, furthermore sensitive to systematic errors.

> Since my samplig-F is not so much higher than the sinewave, I need to
> implement interpolation to get a higher resolution.
> The autocorrelation of a sine wave is a cosine. Steepness around 0&deg; is
> very low. What algo is recommended for interpolation?

Do you know the approximate result to get rid of the many matches of the
sinusoid wave?

> Shall i use a different approach in general?

Use white noise as reference signal rather than sinusoid. More
specifically use white noise with a periodicity of exactly one FFT
length. You can get this easily by taking the inverse FFT transform of
all amplitudes 1 and random phase. This turns the white noise into a
discrete spectrum which might be used to cancel any other, uncorrelated
noise later.

Usually it is also recommended to avoid frequencies close to the Nyquist
frequency of your ADC to avoid artifacts of the anti aliasing filter.
Since at least some parts of the anti aliasing filter are analog
components with tolerances you should be aware that the resulting phase
results at high frequencies. And phase differences of an ordinary 1st
order filter with only one capacitator start far away from the cut off
frequency. This is the second reason to avoid the last ~15% of the
bandwidth. If your ADC is only AC coupled the same applies to the lower
cut off frequency as well. But as the frequency is lower you need not to
be that far away in count of frequency bins.

Simply set the amplitudes of your generated white noise at the unwanted
frequencies to zero 0 rather than 1. The inverse transform will give you
the optimal reference signal for your equipment. Also note that
inaccuracy of the DAC does not impact the result as this is common to
both inputs and will cancel.

To get the final result simply do the following:

Z(f) = FFT(ch1) / FFT(ch2)

I.e.: calculate the FFT of both input channels. Divide the complex
results of each frequency bin. This is the complex impedance of your sample.

Now, if you are interested only in the linear phase shift fit a line
through the phase:

arg(Z(f)) = x * 2 Pi f

You can do this by a chi&#2013266098; method as already mentioned. Fit for the the x
that gives you the least residuals. This is the wanted group delay of
your sample, i.e. the final result.

Of course, you need to discard the frequencies that you have dropped
To be more flexible, you can also use the reciprocal amplitudes at both
inputs as weight factor. Now there is no more need that the reference
signal must satisfy any constraints except for the time synchronicity.
Both methods could be combined, i.e. drop frequencies without energy and
weight the others.

Note: you can get even more output from your measurement, if you use a
more complex model. E.g.:

Zmodel = a * Zmodel1 + b * Zmodel2 + c * Zmodel3

You can still use the chi&#2013266098; method to get a, b and c from your measured
Z(f). From the mathematical point of view it is still a linear fit and
you can calculate it by a pseudo inverse matrix. The method is called
PCA = Principal Component Analysis.

This could either be used to identify different impedance components,
e.g. ESR, ESC and ESL of an electrolytic capacitor, or to separate
systematic errors like parasitic line resistance or capacity.

There are still two systematic errors left in your measurement: cross
talk and channel differences. Almost any ADC has some sort of cross talk
and tolerances, usually due to the analog circuit board. This will
degrade your results, at least at the bleeding edge. But is is easy to
compensate for that.

You need 3 calibration measurements:
#1: ch1 receives the input, ch2 is grounded.
This gives you the ch1 -> ch2 cross talk XT12.
#2: the other way around.
This gives you the ch2 -> ch1 cross talk XT21.
#3: both channels receive the same input.
This gives you the differences between channels D21.

The calculation of the XT12 etc. is the same as above, i.e.
FFT(XT21) = FFT(ch2)/FFT(ch1) ...
At the end you can get a 2 by 2 matrix from this calibration where each
element consists of a frequency dependent table for each FFT bin - in
fact this is a tensor.

| ch1 -> ch1   ch2 -> ch1 |
|                         |
| ch1 -> ch2   ch2 -> ch2 |

There is one thing missing: you have only 3 calibration measurements for
4 matrix elements. As you do not know the absolute input signal you
cannot correct for the absolute response of your inputs. But this is
needless since the absolute input will cancel at the division of both
FFTs anyway. So just set the ch1 -> ch1 factor to all 1.

Now your parasitic transformation matrix looks like:

|     1      FFT(XT21) |
M = |                      |
| FFT(XT12)   FFT(D21) |

i.e. you measure ch1', ch2' instead of ch1, ch2:

| ch1' |   |     1      FFT(XT21) |   | ch1 |
|      | = |                      | * |     |
| ch2' |   | FFT(XT12)   FFT(D21) |   | ch2 |

To get your original values back just take the inverse of M.
Uh, M is a tensor. Don't worry, since time (and frequency) is an
invariant in a linear system you just need to take the inverse at each
frequency bin separately, and well, 2*2 matrix inversion is straight
forward even with complex numbers. And since the determinant cancels in
the result you can just ignore it.

Now you get the corrected results from your measurement:

| FFT(ch1) |                    |  FFT(D21)   -FFT(XT12) |   | FFT(ch1') |
|          | is proportional to |                        | * |           |
| FFT(ch2) |                    | -FFT(XT21)       1     |   | FFT(ch2') |

Only the last division of FFT(ch1) / FFT(ch2) is missing from this formula.
This gives you *very* accurate results even with cheap ADC devices.

I think for only linear phase shifts only the D21 coefficient is
important. But it is straight forward to use the full correction anyway.

Further hints:

- Do not use any window functions. They are not needed because of the

- Strictly avoid overdriving, of course. It would add significant noise
to the measurement.

- Due to the linearity of the process you may add consecutive packages
of samples with the FFT length before further processing. This reduces
noise and also reduces the required CPU capacity because less results
have to be calculated per time at a continuous measurement.
But don't care too much. Even a Raspberry Pi 1 has by far enough power
to do this in real time as long as you run the FFTs on the GPU.

- And do not worry about the number of bits of your ADC. Even an 8 bit
ADC will perform surprisingly reasonable. Using 24 or more bits will

- Last but not least the FFT size. In general it just restricts the
maximum delay you can safely measure. So taking a larger FFT size will
not improve results in general. A smaller FFT is just equivalent to
dropping each frequency except for every n-th one in the reference signal.
There is no direct dependency to accuracy as long as you are talking
about linear phase shifts only. But you become more sensitive to noise
and/or systematic errors at certain frequencies. So don't go to far. I
would recommend 256 at least.

- If you are in worry about small non linearities in your system just
drop every even frequency bin from your reference signal. This will
result in cancellation of any second order harmonic in the result
because their energy will be at sums or differences of two other
frequencies each. This sums or differences will be necessarily at /even/
frequency bins which will be ignored. The drawback is that it will
double the measurement time to get the same result accuracy.
Unfortunately it is not that easy to cancel higher order harmonics, e.g.
H3, from the result. At least I did not find out how to do this so far -
except for measuring each frequency alone which takes ages, of course.

Marcel

On Monday, July 31, 2017 at 1:23:09 PM UTC-5, herrman...@gmail.com wrote:

> One possibility is to do a non-linear least squares fit to
> A*sin(B*t+C)  If B is well enough known, and since you are
> generating the sine it should be, you could do a linear least
> squares fit to D*sin(w*t) + E*cos(w*t).  (That is, linear in
> the coefficients D and E.)  You can do this even if the sampling
> rate is very low, if you have enough samples.

"EVALUATION OF A STANDARDIZED SINE WAVE FIT ALGORITHM", Peter Handel, Royal Institute of Technology, Stockholm.

"Standardized sinus-wave fitting algorithms, extensions and applications", Tomas Andersson and Peter Handel, Royal Institute of Technology, Stockholm.

Thank you Marcel & Greg.
I think I've got enough material for this project and two generations of
improvement. :).
I will come back on you all when I've implemented some of the ideas in
this topic.
Best regards.
---------------------------------------
Posted through http://www.DSPRelated.com

On Monday, July 31, 2017 at 9:38:26 PM UTC-7, robert bristow-johnson wrote:

(snip)

> your sample rate still needs to be faster than twice that of the sinusoid.

No.

This is almost the perfect case for bandpass sampling.
(That is, with a bandwidth of zero.)

If you have a sine of known frequency, but unknown phase and
amplitude, two samples are in theory enough.  If things aren't
quite that perfect, you should have more samples.  For one,
you can average out noise, which might some in.

Next, consider the case where you know the fundamental
exactly, but not the spectrum. You could sample a signal
of period T every 10.001 T for thousands of samples.
With just 1000 samples, you have the equivalent of 1000
samples over one period.

(snip)


benjamin.couillard@gmail.com writes:

> If you only need to compare the phase, you could use the Hilbert
> transform on each of the signals. Extract the instantaneous phase from
> each signal. Subtract the phase and you will get the phase difference
> between the 2.

I don't know if it's optimal, but I sure like this approach
conceptually.