I'm not sure if this is the best forum for this or whether stackoverflow mightn't be better, but I'll ask here and see what people think.

I have a system measuring 50Hz voltage and currents at 1kHz. I'm using quadrature demodulation to measure the phase shift and calculating cos(φ) from that, and also calculating the active and RMS power and inferring cos(φ) from them. These calculations are done over one-second slices of the data.

I'll attach example Python code below but basically the method looks like this:

  • Multiply the instantaneous voltage and current timeseries together to get instantaneous power, p_inst.
  • Integrate p_inst (using a trapezoidal integration) and divide by the time length of the data slice to get the active power, p_act.
  • Square each sample in p_inst, take the average of p_inst**2 across the timeseries and then take the square root of that value to get the RMS power, p_rms.
  • Calculate the power factor, pf = p_act / p_rms.
  • Using a zero-crossing-detection method, estimate the frequency of the voltage signal.
  • Generate quadrature reference signals at that frequency (ie cos(2πft), sin(2πft)).
  • Multiply the voltage by each of the reference signals and integrate (again a trapezoidal approximation) each of the resulting timeseries, resulting in v_cos and v_sin. The phase angle of the voltage signal, v_phi, relative to some arbitrary reference, is then atan2(v_sin, v_cos).
  • Repeat the previous step for the current. The phase lag of the current relative to the voltage is then φ = v_phi - i_phi.
  • My second estimate of the power factor is then cos(φ).

For a sinusoidal signal with 30 degrees phase lag, the correct cos(φ) is approximately 0.86602540378443871. The quadrature demodulation method produces a very good approximation, 0.86636025346085943. But the ratio-of-powers method produces a very wrong estimate - 0.77542956418409648. This is equivalent to nearly ten degrees error in the phase angle.

I assumed at first that I'd got the quadrature demodulation wrong (being the more complex calculation) but it produces the right answer. I then assumed that the signal was badly non-sinusoidal and that this would explain the difference - but the code below does the same calculation on ideal sinusoids.

What have I got wrong here?

Full Python code demonstrating the problem:

import numpy as np, pandas as pd
from matplotlib.pyplot import *

t = np.arange(0, 600, 0.001)
t_rad = 2 * np.pi * 50 * t
phi = 30 * np.pi / 180

v = 230 * np.sin(t_rad)
i = 200 * np.sin(t_rad + phi)

p_inst = v * i

data = pd.DataFrame(np.array([t, v, i, p_inst]).transpose(), columns=['t', 'v', 'i', 'p_inst'], index=t)

def gen_act(x):
    return np.trapz(x.p_inst, x=x.index) / (x.index[-1] - x.index[0])

def gen_rms(x):
    return np.sqrt(np.mean(x.p_inst**2))

second_bins = data.groupby(lambda x: int(x))
p_act = second_bins.apply(gen_act)
p_rms = second_bins.apply(gen_rms)

def estimate_frequency(t, v):
    t = t.values
    v = v.values
        zero_crossings = np.where(np.diff(np.sign(v)) > 0.5)[0]
        diffs = np.diff(t[zero_crossings])
        accum_intervals = np.copy(diffs)
        accum_interval = 0
        ii = 0
        for ii in range(len(accum_intervals)):
            accum_interval += diffs[ii]
            accum_intervals[ii] = accum_interval
            if accum_interval >= 0.01:
                accum_interval = 0

        zero_crossings = zero_crossings[np.hstack([np.where(accum_intervals > 0.01)[0], -1])]
        ends = zero_crossings[[0, -1]]
        frequency = (len(zero_crossings) - 1) / (t[ends[1]] - t[ends[0]])
        return frequency
        return float('NaN')

def generate_reference_signals(t, f):
    freq = np.mean(f)
    if np.isnan(freq):
        return None
    x = t * 2 * np.pi * freq
    return np.sin(x), np.cos(x)

def estimate_phi(t, v, i, refs):
    def angle_from_refs(t, x, refs):
        s_i = np.trapz(x * refs[0], t)
        c_i = np.trapz(x * refs[1], t)
        return np.arctan2(c_i, s_i)

    angle = angle_from_refs(t, i, refs) - angle_from_refs(t, v, refs)
    while angle > np.pi:
        angle -= 2*np.pi
    return angle

def gen_phi(x):
    f = estimate_frequency(x.t, x.v)
    refs = generate_reference_signals(x.t, f)
    return estimate_phi(x.t, x.v, x.i, refs)

phi = second_bins.apply(gen_phi)
cos_phi = np.cos(phi)
pf = p_act / p_rms

print('Power factor estimated from quadrature demodulation: {}'.format(np.mean(cos_phi)))
print('Power factor estimated from measured power: {}'.format(np.mean(pf)))
print('True power factor: {}'.format(np.cos(30 * np.pi / 180)))

2 Answers 2


Your third step, squaring p_inst, is incorrect.

What you really want is to find the RMS values of the sets of voltage and current samples separately, yielding v_rms and i_rms, and then multiply those results together to get p_rms.

  • \$\begingroup\$ Ergh, aren't they the same thing? \$\endgroup\$
    – Tom
    Jul 22, 2017 at 11:39
  • \$\begingroup\$ No, they are not, and that's the whole point. \$\endgroup\$
    – Dave Tweed
    Jul 22, 2017 at 12:35
  • \$\begingroup\$ In my defence, the bug is one I inherited. But I should have spotted it. \$\endgroup\$
    – Tom
    Jul 22, 2017 at 12:39

This is basically a naming issue, which once more proves how important is using correct vocabulary in technical matters.

You say you are after RMS power, then your original algorithm is ok since given N instantaneous power samples


the above is just applying the definition of Root of Mean of Squares.

Now the point is RMS power is a mostly useless quantity and definitevely not what you need to work out power factor.

What you should have computed is apparent power instead, often designed by capital S, and [S]=VA.

One possible generalisation of apparent power which works well even with non sinusoidal distorted currents and voltages simply states that


and also, calling Pact the average, i.e. active power $$P_\text{act}=\text{PF}\;S\quad\Rightarrow\quad \text{PF}=\frac{P_\text{act}}{S}=\frac{P_\text{act}}{V_\text{RMS}I_\text{RMS}}$$

Note that power factor \$\text{PF}=\cos\varphi\$ in case of sinusoidal quantities but also includes non linear effects in case there are some.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.