Calculating xPower

[Steve Mansfield]

G'day

I'm trying to replicate the xPower calculation, but not having any luck.

Not being familiar with exponential moving averages, I did the usual Google search, and found this page with a sample: http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:moving_averages, and they have a sample file to try:  http://stockcharts.com/school/data/media/chart_school/technical_indicators_and_overlays/moving_averages/cs-movavg.xls

If I understand correctly, the time period for xPower is 30s, so the smoothing factor is 2/31 = 0.06451612903, and so the next EMA would be

=(2/31)*(B3-C2)+C2

if power is in column B, and I want the EMA to be in column C.

So I loaded a fairly small file into GC, and copied the Power column into a Sheet, and applied the above formula, then taking each EMA , raising to the 4th power, average of all those, and fourth root.

Thing is that I get 225W, GC gets 215, so I'm doing something wrong

Here's the Sheet if you fancy trying yourselves: https://docs.google.com/spreadsheets/d/1BdgCjlycvYIf_3KMObOLK2BcrUfQV-RAZf8pF4hdDL0/edit?usp=sharing

Any ideas?

Thanks!

Steve

[Mark Liversedge]

Here is the code for the xPower metric

       static const double EPSILON = 0.1;

        static const double NEGLIGIBLE = 0.1;

        double secsDelta = item->ride()->recIntSecs();

        double sampsPerWindow = 25.0 / secsDelta;

        double attenuation = sampsPerWindow / (sampsPerWindow + secsDelta);

        double sampleWeight = secsDelta / (sampsPerWindow + secsDelta);

        double lastSecs = 0.0;

        double weighted = 0.0;

        double total = 0.0;

        int count = 0;

        RideFileIterator it(item->ride(), spec);

        while (it.hasNext()) {

            struct RideFilePoint *point = it.next();

            while ((weighted > NEGLIGIBLE)

                   && (point->secs > lastSecs + secsDelta + EPSILON)) {

                weighted *= attenuation;

                lastSecs += secsDelta;

                total += pow(weighted, 4.0);

                count++;

            }

            weighted *= attenuation;

            weighted += sampleWeight * point->watts;

            lastSecs = point->secs;

            total += pow(weighted, 4.0);

            count++;

        }

        xpower = count ? pow(total / count, 0.25) : 0.0;

        secs = count * secsDelta;

        setValue(xpower);

        setCount(secs);

[Steve Mansfield]

Thanks, so firstly I was using NP's 30s not 25s per Skiba.

If I do that, the xPower result becomes 227W which is still different.

[andy aardema]

In your spreadsheet, here's the formula I see for weighted average:

=((2/C$1)*(B3-C2))+C2

I'm not sure where the "2" above is coming from.  Each new value should be the current power * (1 / 26) plus the previous weighted value * (25 / 26).  (In Mark's code, for 1 second samples, sampsPerWindow should evaluate to 25, attenuation will be 25 / (25 + 1), and sampleWeight will be 1 / (25 + 1) ).

I've implemented xP in this way writing a connect iq app for my garmin edge head unit, and the overall xP values for rides match exactly to GC when I import.

I made a copy of your google sheet and changed the weighted average column to this:

=(1/26*B3)+(C2*25/26)

That yields an xP value of 215.4.  The other issue is the very first entry for weighted average shouldn't be =B2.  It needs to be =(1/26*B2)+(0*25/26).

Andy

[Steve Mansfield]

Thanks!

That all makes sense. The "2" just comes from the example I found.

From https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average it looks like we are using a modified moving average where  ?

Here's where I get confused. From the BikeScore paper figs1 and 2:

In fig1, the EWMA (red) tracks more closely with instantaneous power than SMA (blue), which is the purpose of EWMA. In fig2, EWMA is slower to respond than SMA.

From Wikipedia: "An exponential moving average (EMA), also known as an exponentially weighted moving average (EWMA),^[5] is a type of infinite impulse response filter that applies weighting factors which decrease exponentially. The weighting for each older datum decreases exponentially, never reaching zero"

Here's my example ("25" is the formula I got from the Stock Market example), "EWMA" is the aforementioned 25/26 algorithm, and SMA is a straight 25s moving average.

I confess to being confused about Skiba's figs 1 and 2. Could it be that in fig2 the colours for EWMA and SMA are reversed?

[Steve Mansfield]

Quick note: in my chart the horizontal axis is called "Pwr", its not, it the datum number.

[andy aardema]

I think Skiba's plots are right.  One issue is his x-axis is much more condensed.  It might look like the simple moving average "tracks" instantaneous power better, but the idea is it doesn't track real exertion as well.  If you think about your heart rate during a sustained, maximum-intensity interval followed by complete coasting, it trails off very much like the red line in Skiba's fig. 2.  A simple weighted average would imply that your heart rate would drop right to resting after 30 seconds, when it actually decays over many minutes (i.e. exponentially, and in theory out to infinity).

Here's a screenshot that I find pretty neat.  It's a plot of my commute to work (zero power is at stoplights).  The agreement/alignment between xP and HR over the short term is rather amazing.

Andy

[Steve Mansfield]

Interesting, here's a similar analysis from a ride earlier today.

To get "normalised" I've used the more usual mathematical definition ie n/average(n). Because HR doesn't go to zero I used (n-min(n))/(average(n)-min(n)) for HR.

[Steve Mansfield]

OK, final word from me on this, since I can now replicate the xPower calculation.

If the aim of using an exponential moving average with a lag of 25s is to model oxygen kinetics etc, the formula used ((1/26*B3)+(C2*25/26)) sure does lag a lot. The more regular means of calculating EWMA mimics a 25s response much better.To illustrate, here's the first 10 minutes of a road race. I've use "EWMA" to mean the more usual formula I mentioned above.

Interestingly in Skiba's paper the EWMA moves more closely in the one diagram, and less closely in the other.

[andy aardema]

Heh, I'm still following you and still interested; I'm just busy with other things.  ;-)

Yes, I think the lagginess is actually a feature and/or by design.  It trails actual power by a good margin, but it tracks heart rate.  At least that's what I typically see in my plots.  The Y axis would need some pretty fancy scaling to mate HR and power, but along the X axis the vertical peaks and valleys line up extremely well for me.  They're usually within a second or two during highly variable riding, which I take to mean the decay time-constants seem about right (and xP trends to my HR slope better than the linear decay of NP).  That's just one piece to tracking exertion/training response, and I would add that I'm not sure that the many fancy constructs really do any better than avg. power * time over the big picture.  From a macro view, they all look pretty similar.  I think Mark L. has made that point before and posted graphs to illustrate it.

Related to traditional EWMA calculation, that has me curious.  

I thought traditional was along the lines of:

 s(t) = alpha * p(1) + (1 - alpha) * s(t - 1)

which is what we're doing.  

My framework for EWMA is signal processing, where it's a standard infinite impulse response filter.  In the case of IIRs, they by definition have a phase delay seen as the timing difference between sinusoidal in => filter = sinusoidal out, and that's the lag on your (very cool!) graph.

Andy

[Steve Mansfield]

I think we agree on the principle of using EWMA to model physiological response, and I'm pretty sure Skiba knows more than me :-)

So the only question I have left is if we're using the correct lag / smoothing factor.

Reason I ask is that when I searched for how to calculate EWMA, those sites (usually stock market / financial market related) such as this one, that have worked examples where K = the smoothing factor, define as

k=2/(D+1)

where D = the number of days, or for us the number of seconds.

The formula we've established Skiba uses, and replicated by GC is

k=1/(D+1)

How does this change things, other than make xP lag more? The above website would say that in essence Skiba's algorithm uses a smoothing period of 51s, rather than 25s.

For fun, I imagined a session of intervals, where the cyclist produces 1000W for 120s, rests for 120s and so on. Because xP uses k=1/26 it means it takes nearly 2 minutes for xP to almost catchup.

[andy aardema]

I had a bit more time to study those interval graphs.  The first thing I'd note is that, while your gold xP line rises more slowly, it falls more slowly as well.  Average may be a wash.   

Next, your plot lines over 120 seconds intervals were not quite what I expected.

I mocked up an excel graph with the same 2 min power intervals at 1000 watts.  I've done the full avg( {set}^4 )^0.25 operation for NP and xP, and I added 30 second windows for both as well as average for ride duration.

These values match what I get in GC and my own garmin app, but I was left puzzled on the climb/decay rate of xP.  The "half-life" doesn't appear to be 25-26 seconds (or 51), before or after the 4th order polynomial operation.  It's closer to 18.

I dug through formulas for exponentially smoothed averages versus an actual exponential decay:

https://en.wikipedia.org/wiki/Exponential_smoothing

https://en.wikipedia.org/wiki/Exponential_decay

The EWMA formula I've been using is s(t) = alpha * x(t) + (1 - alpha) * s(t-1)

where alpha = 1 / 26.

I think we've all thought that meant the "half-life" was roughly 26 seconds (at least I did).  According to the first link, alpha is actually the smoothing factor.  I really wanted to know the relationship between lambda and tau of traditional exponential decay equations and alpha, so down the rabbit hole I went.  It turns out it's not that hard.

Again from the first link, exponential smoothing is typically discrete and decay continuous, so just to be precise: 

alpha = 1 - e ^ (-1 / tau)

where tau is the mean lifetime = time constant, and should in practice by very close to alpha when the samples delta T is small compared to alpha.  Some rearranging there, and I have

-1 / tau = ln(1 - alpha)

or

tau = 25.49673, so that's pretty close.

Half-life is simply tau * ln(2), or 17.67 seconds.

Going back to my spreadsheet, indeed my "25 second EWMA" column reaches 500 watts after 18 seconds of the 1000 watt interval, and likewise drops by half 17-18 seconds after an interval ends.  The actual xP column is more complicated due to the 4th order averaging.


As to whether that means anything, I don't think so.  I believe the original Skiba paper indicates that they picked a time constant which matched the decay plot of actual measured physiological markers and not the other way around.  In any case, it's an interesting exercise.  Separately, one can see how NP ends up higher than xP over equivalent intervals.  That's always the case in my ride files, and frequently with hard variable rides I find maximal 20 min NP values impossibly high (i.e. well above my mean maximal power curve), while maximal 20 min xP values tend to be right at the limit.  That lead me to trust xP more, but perhaps it's just been luck.  The EWMA of values I get. The (^4)^0.25 part, as I understand it, attempts to weight variability higher and/or account for the non-linear costs of working well above your threshold for sustained durations.  That also makes intuitive sense, but how can it not by anchored by one's own FTP or CP (or mean maximal curve?)?  Riding intervals alternating from, say, 0.7 FTP on and 0.1 FTP off means nothing compared to doing the same steady-state average.  Riding intervals at 1.4 FTP on and 0.8 FTP off means a lot.  Those two "variabilities" are the same under this model, when in practice variability means a whole lot more at/above threshold than it does below it.  That's my interpretation, anyway.  Perhaps I'm missing something fundamental.

Andy


On Thursday, September 15, 2016 at 4:40:07 PM UTC-4, Steve Mansfield wrote:

[Steve Mansfield]

This has all been interesting... not sure where it leaves us!

Here's a chart from a recent set of intervals where I was attempting to hold 400W for 6 minutes. I've not plotted the values raised to the 4th power, averaged etc because the paper says that doing so for anything less than 20 minutes gives unreliable results. I've scaled HR on the r/hand axis try to follow the power curves.

As before EWMA = the formula I found on finance web sites, xP is Skiba's, nP is Coggan's

[mickebergma...@gmail.com]

After reading these posts I added xPower to the metrics I use. And from my side it correlates with feeling way better than NP.
I.e. yesterday I did an interval training on the Wattbike with AP 207 W, NP 264 and XP 236.
As 264 is over my FTP it reflects really badly the feeling from the (relatively easy) session. XP feels OK.

Then a question. I know I could search the FAQ, Wiki, Etc but I'm not nerd enough to find it even if I spent some time trying.

Is TriScore based on NP or XP? It is about the same amount lower than TSS as XP from NP (and hence, also reflects feeling better than TSS).
If not what's the difference in the calculation.