Forums

Sine Wave autocorrelation, interpolation of phase

Started by eduardoG26 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 µ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° 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 µ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° 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 µ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° 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 µ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 differences of your two inputs will directly impact your measurement 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 from your reference signal. 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 synchronous measurement and will significantly degrade your results. - 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 probably mainly add meaningless noise. - 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. Additionally, note the following: 1. You will have to wait some time (depending on the filters used for the Hilbert transform) for the startup transients to die out. 2. If you have N post-transient samples available, you can (and should) average them to improve the phase estimate. -- Randy Yates, DSP/Embedded Firmware Developer Digital Signal Labs http://www.digitalsignallabs.com