Phase angle calculation

[kvresto]

Hi Everyone.

I’ve been trying to calculate the phase angle between two sine waves differing only in phase, and I’ve got it working in MATLAB, using the FFT of each sine wave and then applying MATLAB’s angle function. (as an example I think from the mathworks site??)

 Now I’m trying to do  the same with the STM32F3. I’m using the DSP Lib that can be downloaded from ST, and I’ve gotten as far as the FFT part, but as I’m not at all familiar with the library, I’m hoping someone with more experience can help with some direction as to how to get the answer I’m looking  for. Unfortunately I'm also a bit hazy on the math's, so not starting from a strong base, but the Lib should be able to do most of the hard work.

Any ideas??

thanks
kvresto

[coppice]

You might want to seek help with a much simpler and more accurate approach before you look for help with your current solution.

One simple and effective approach is to take the inner product of the two waveforms, and the inner product of the two waveforms with one delayed by a well known amount (e.g. a whole number of samples). You now have two measures of signal correlation. If the two waveforms are exactly in phase you will see the maximum possible correlation. If they are exactly in quadrature you will get zero correlation. Apply a little trigonometry, you will get your answer with good precision.

[JPortici]

F3 has a cortex M4 with FPU, right? if it's the case you ABSOLUTELY have to use it. 512 pt fft done on a 24 MHz cortex m3 took like 100 ms

Anyway, you should look at the cmsis manual. if i remember correctly the fft routine returns a structure with two series of data: real and imaginary bins from 0 to Fs (or, if you want.. first is DC then for the rest the first half is from dc to Fs/2 and the other half to -FS/2 to dc)
then you find out the phase in the usual manners

Question: you say quote "two sine waves differing only in phase"
i assume that the module and the frequency is identical.. can't you find out the phase by other simpler means? (academic purpose maybe?)

[kvresto]

coppice, thanks for your advice, I will look into the method you described.

JPortici, your comment was:

then you find out the phase in the usual manners

can you go into a bit more detail on how to do this academically, and using the DSP Lib?

Also which FFT function are you referring to, is it  arm_rfft_fast_f32()?

[JPortici]

that one, if you want to work with floating point (and you want to)

the result will be of two arrays in the form a[n] + i*b[n], where n is the frequency bin. from that, finding out the phase is trivial (you should also have ready made functions in the cmsis as well. again, refer to the cmsis documentation)

[KE5FX]

Lots of ways to skin this particular cat.  Rather than a full-blown FFT, you might consider running the sines through Hilbert transforms, and  taking atan2(q,i) on the resulting analytic signals.

[DanielS]

I like the dual convolution method: if the two signals are the same amplitude and frequency (simply normalize both to eliminate amplitude as a concern) then the result will be 0.5 if the two waveforms are in-phase, -0.5 if they are 180° out of phase, 0 if they are +/- 90° out of phase, sqrt(1/8) for +/- 45°, etc. You need the second convolution with a delayed waveform to discriminate which quadrant you are working in.

[kvresto]

thanks all for your comments.

DanielS the waveforms will not be the same freq and ampl will change as per the application.

I'm using matlabs angle function, which is actually the atan2() function. How would this same process be done in my case(DSP_Lib), and please don't say its "trivial" again


PS: if your not sure, just say so, were all just learning here, I know I am.

[coppice]

I like the dual convolution method: if the two signals are the same amplitude and frequency (simply normalize both to eliminate amplitude as a concern) then the result will be 0.5 if the two waveforms are in-phase, -0.5 if they are 180° out of phase, 0 if they are +/- 90° out of phase, sqrt(1/8) for +/- 45°, etc. You need the second convolution with a delayed waveform to discriminate which quadrant you are working in.

I assume you mean dual correlation, not dual convolutuion, and are referring to what I wrote earlier.  What you describe would give very poor accuracy for a wide range of phase shifts. First, two signals would rarely be the same amplitude after some manipulation. A practical solution would need to work out their relative gains with, say, a pair of RMS measurements conducted in parallel with the correlations. You would have very poor accuracy when the signals are nearly in phase, as the difference in correlation between signals in phase and signals one or two degrees out of phase is very small. As the phase difference increases the measurement accuracy would improve. Above 20 degrees or so, the measurement would be pretty good. The whole point of using delayed and undelayed signals is to avoid this. With the results of the two correlations, and some trigonometry, you will get a fairly balanced level of accuracy for any phase shift, and the relative gains of the two channels doesn't matter.

[kvresto]

Also KE5FX, how is what you have suggested too different from a straight up FFT?, its the atan2 part I'm looking to solve using this Lib for anyone with experience, ie how to format my fft results so as to plug them into atan2?    Its probably trivial!

[KE5FX]

Also KE5FX, how is what you have suggested too different from a straight up FFT?, its the atan2 part I'm looking to solve using this Lib for anyone with experience, ie how to format my fft results so as to plug them into atan2?    Its probably trivial!

The Hilbert transform would likely be implemented with an FIR kernel, which might or might not outperform an FFT of a given size depending on how the pipelining is done.  It's just two ways of doing the same thing.

I'd suggest having a look at this thread on comp.dsp, as well as Googling "polar discriminator."

[coppice]

Also KE5FX, how is what you have suggested too different from a straight up FFT?, its the atan2 part I'm looking to solve using this Lib for anyone with experience, ie how to format my fft results so as to plug them into atan2?    Its probably trivial!

The Hilbert transform would likely be implemented with an FIR kernel, which might or might not outperform an FFT of a given size depending on how the pipelining is done.  It's just two ways of doing the same thing.

I'd suggest having a look at this thread on comp.dsp, as well as Googling "polar discriminator."

The effectiveness of a Hilbert transform, implemented by the traditional FIR approach, very much depends on where the frequency of the signal of interest sits in the band. Just a few FIR terms is enough to build a pretty accurate transform if the signal is near the middle of the band, but it takes a huge number of terms if the signal is near one end of the band. You can implement a pseudo-Hilbert using 2 IIR filters, which shift the signal in two directions, so their mutual outputs are 90 degrees apart at all frequencies. So, all is not lost if the signal is near one end of the band, but you need to take considerable care. You will still need to perform correlation between the analytic signals from the Hilbert transforms, and do trigonometry with the results, or the answer will be very susceptible to noise.

[wolf32d]

Hi Everyone.

I’ve been trying to calculate the phase angle between two sine waves differing only in phase, and I’ve got it working in MATLAB, using the FFT of each sine wave and then applying MATLAB’s angle function. (as an example I think from the mathworks site??)

 Now I’m trying to do  the same with the STM32F3. I’m using the DSP Lib that can be downloaded from ST, and I’ve gotten as far as the FFT part, but as I’m not at all familiar with the library, I’m hoping someone with more experience can help with some direction as to how to get the answer I’m looking  for. Unfortunately I'm also a bit hazy on the math's, so not starting from a strong base, but the Lib should be able to do most of the hard work.

Any ideas??

thanks
kvresto

just calculate the xcorr (cross correlation) between the two functions and take the absolute maximum of the result

[kvresto]

wolf32d, thanks, how do you then define the lag in terms of degrees phase shift?

I don't know what the "absolute maximum " is all about?

Have you tried this with signals of different amplitude, but with the same freq?

[mark03]

DanielS the waveforms will not be the same freq and ampl will change as per the application.

If the two waveforms aren't the same frequency, the phase angle is [obviously] a continuously changing quantity.  How far apart in frequency are we talking here?  Most of the proposed methods assume an average over many cycles...

[kvresto]

mark03, its an (experiment)cct, that will measure ESR, and impedance of a capacitor(for now).

Series RC with voltage measured across Rs and Cs. Freq is fixed, but ampl will change.

cheers

[kvresto]

Also thanks everyone for some very good information. your help is very welcomed.
Matlab is now running, and I've got my old text books out, having fun.

kvresto.

[wolf32d]

wolf32d, thanks, how do you then define the lag in terms of degrees phase shift?

I don't know what the "absolute maximum " is all about?

Have you tried this with signals of different amplitude, but with the same freq?

There you go:

Code: [Select]


% generate some signals
sampling_step = 0.01;

t = 0 : sampling_step : 4*pi; %generate a timebase

samples_per_period = 2*pi / sampling_step;

v1 = sin(t);

phase = pi/5; % manually introduce an arbitrary phase shift

v2 = sin(t-phase);

plot(t, v1, t, v2);

%%%%%%%%%%%%%%%%%%%%%%%% now he have v1 and v2 i.e. our simulated signals

correlation = xcorr(v2,v1); %evaluate cross correlation

[corr_max, corr_argmax] = max(correlation); %evaluate the max of correlation and its position

%plot(correlation);

measured_phase_in_samples = corr_argmax - length(correlation)/2; % evaluate distance of max position from the center

measured_phase_angle = (2*pi / samples_per_period) * measured_phase_in_samples % convert to fraction of period (2 pi) i.e. agle in radians
                                                                               % should be equal (or very close to) "phase"

%evaluate phase in degrees
measured_phase_angle_deg = 180 * measured_phase_angle / pi

if you have real signals you should evaluate "samples_per_period" by other means
e.g. if you sample a 10Hz signal at a sampling rate of 1kHz you get samples_per_period = 100

[kvresto]

if you have real signals you should evaluate "samples_per_period" by other means

good example wolf32d. Something like 1/f *Fs should do it. EDIT: better yet decimate over a larger sample range.

Also the result is about 1.5 degrees off, seems close enough, can it be improved? (I'm still testing with it).

kvresto

[coppice]

The FFT, tranforms and correlation based solutions are like using a big bertha canon to kill a mosquito. And I'm choosing my words very carefully.

Your signal 1 with phase alpha1 is y1 = sin(x-alpha1), and the signal 2 is y2 = sin(x-alpha2).

Let's use some trigonometry! We take the product : = sin(x-alpha1)*sin(x-alpha2)

We can transform this into a sum using the trigonometric identities : = 0.5 * ( cos(alpha1-alpha2) - cos(2*x - alpha1 -alpha2) )

If you take the mean over enough periods, the term in cos(2*x... becomes negligible, and you're left only with the term in cos(delta_phase). It's possible to calculate over a low count a periods too : If you know the signal frequency beforehand, you can take the mean over an integer number of periods.

A Matlab code is attached. For two samples of length "n", there's only n multiplications.

Your approach assumes the gains of the two signals are identical, which is unlikely to be the case. Its also assumes you aren't looking for a great deal of accuracy, because the accuracy is going to vary a lot with the size of the phase shift. The correlation approach only requires 2N multiples, and solves both issues.

[kvresto]

EDIT:

[wolf32d]

On my computer, computing 1000 times my example code takes  0.060s, and the correlation approach (Exact code given by wolf32d, giving in argument the same vectors) takes 1.975s.

I definitely agree for the amplitude dependance. But assuming A1 is the amplitude of signal 1 and A2 the same for signal 2, you're left with a term in A1.A2.cos(delta_phase). Depending on your application it still may be relevant to measure both of the amplitudes and divide.

LOL WHAT WAS I THINKING? XCORR SCALES o(len^2)

I'm very attached to the xcorr function (I used it a lot to calculate delays in complicated signals where only correlation helps)... but your solution is actually VERY GOOD for sinusoidal stuff,
besides, if the two signals have different gains and offsets this transform will do the trick  (works a treat also with noisy stuff):

Code: [Select]

% REMOVE OFFSET
y1 = y1-mean(y1)
y2 = y2-mean(y2)

% divide by rms and sqrt(2) i.e. NORMALIZE TO 1 (0 to peak)
y1 = y1/(sqrt(2) * std(y1))
y2 = y2/(sqrt(2) * std(y2))

conclusion:

trigonometric method:
-pro: fast
-con: a bit unreliable if signals are not stable in terms of gain and offset or are nonlinearly distorted

xcorr method:
-pro: robust if signals have harmonic distortion & shit
-con: very slow