Amplitude reduction of Simple Moving Average on sinusoid

Question

When a -causal/non-centered- simple moving average filter (arithmetic mean) with length $n$, is applied to a sinusoidal curve with period $p$, and the resulting curve is compared to the original input signal (the sinusoid), 2 things can be observed:

1. The phase is shifted by $360*((n-1)/2)/p$ degrees
2. The amplitude is reduced.

How do you calculate the amplitude reduction?

Please give a numerical example: say $n = 10$, $p=40$, the amplitude of the sinusoid=$1$. My measured result is $0.8984644$* (the maximum of the red curve on the image below), but can you compute this analytically?

Thank you.

---

*It was actually 0.90124: please see below to the kind answers why this was wrong: Dan Boschen: "the output simply wasn't sampled right at its max value".

Answer 1

Okay, this takes a bit of algebra, Euler's formula, and the geometric series summation formula, and some plugging and chugging, but here is how you can calculate it directly:

$$ \begin{aligned} x[m] &= \frac{1}{n}\sum_{k=0}^{n-1} A \cos \left( (m-k) \frac{2\pi}{p} + \phi \right) \\ &= \frac{1}{n}\sum_{k=0}^{n-1} A \left[ \frac{e^{i\left( (m-k) \frac{2\pi}{p} + \phi \right)} + e^{-i\left( (m-k) \frac{2\pi}{p} + \phi \right)} }{2} \right] \\ &= \frac{A}{2n}\left[ \sum_{k=0}^{n-1} e^{i\left( (m-k) \frac{2\pi}{p} + \phi \right)} + \sum_{k=0}^{n-1} e^{-i\left( (m-k) \frac{2\pi}{p} + \phi \right)} \right] \\ &= \frac{A}{2n}\left[ \left( e^{i\left( m \frac{2\pi}{p} + \phi \right)} \right) \sum_{k=0}^{n-1} e^{i\left( -k \frac{2\pi}{p} \right)} + \left( e^{-i\left( m \frac{2\pi}{p} + \phi \right)} \right) \sum_{k=0}^{n-1} e^{-i\left( -k \frac{2\pi}{p} \right)} \right] \\ &= \frac{A}{2n}\left[ \left( e^{i\left( m \frac{2\pi}{p} + \phi \right)} \right) \sum_{k=0}^{n-1} e^{ \left( -i \frac{2\pi}{p} \right)k} + \left( e^{-i\left( m \frac{2\pi}{p} + \phi \right)} \right) \sum_{k=0}^{n-1} e^{\left( i \frac{2\pi}{p} \right)k} \right] \\ &= \frac{A}{2n}\left[ \left( e^{i\left( m \frac{2\pi}{p} + \phi \right)} \right) \left( \frac{1 - e^{ \left( -i \frac{2\pi}{p} \right)n} }{1 - e^{ \left( -i \frac{2\pi}{p} \right)} } \right) + \left( e^{-i\left( m \frac{2\pi}{p} + \phi \right)} \right) \left( \frac{1 - e^{ \left( i \frac{2\pi}{p} \right)n} }{1 - e^{ \left( i \frac{2\pi}{p} \right)} } \right) \right] \\ &= \frac{A}{2n} \left( \frac{e^{ \left( i \frac{\pi}{p} \right)n} - e^{ \left( -i \frac{\pi}{p} \right)n} }{e^{ \left( i \frac{\pi}{p} \right)} - e^{ \left( -i \frac{\pi}{p} \right)} } \right) \left[ \left( e^{i\left( m \frac{2\pi}{p} + \phi \right)} \right) \left( \frac{ e^{ \left( -i \frac{\pi}{p} \right)n} }{ e^{ \left( -i \frac{\pi}{p} \right)} } \right) + \left( e^{-i\left( m \frac{2\pi}{p} + \phi \right)} \right) \left( \frac{ e^{ \left( i \frac{\pi}{p} \right)n} }{ e^{ \left( i \frac{\pi}{p} \right)} } \right) \right] \\ &=\frac{A}{2n} \left(\frac{\frac{\sin\left(\frac{\pi}{p}n\right)}{2i}}{\frac{\sin\left(\frac{\pi}{p}\right)}{2i}}\right) \left[ e^{i\left( m \frac{2\pi}{p} + \phi - \frac{\pi}{p}( n - 1 ) \right)} + e^{-i\left( m \frac{2\pi}{p} + \phi - \frac{\pi}{p}( n - 1 ) \right)} \right] \\ &=A \left(\frac{\sin\left(\frac{\pi}{p}n\right)}{n\sin\left(\frac{\pi}{p}\right)}\right) \cos \left( m \frac{2\pi}{p} + \phi - \frac{\pi}{p}( n - 1 ) \right) \\ \end{aligned} $$

You can see that this is the original signal equation with a phase adjustment and an amplitude adjustment.

Plugging in your values:

$$ n = 10, p = 40 $$

Into the amplitude adjustment term, you get:

$$ \frac{\sin\left(\frac{\pi}{p}n\right)}{n\sin\left(\frac{\pi}{p}\right)} = \frac{\sin\left(\frac{\pi}{4}\right)}{10\sin\left(\frac{\pi}{40}\right)} \approx \frac{0.7071}{10 \cdot 0.07846} \approx 0.90124 $$

Which differs from your answer somewhat. The mistake could be mine, I did this quickly.

Side note to Hilmar and others: The sinc function does not apply here. It is only an approximation in this situation. You need to use the discrete sinc function which is also called the alias sinc function or the Dirichlet kernel.

---

Followup from Ed V's answer:

As I thought, the discrepency lies in the sample points not falling at the peak.

Ed V's answer clearly shows the continuous case is the limit of the discrete case as the sampling density increases. His averaged functions aren't the same though.

I wrote a quick program to find the peak point and approximate the averaged signal with a parabola. The peak ends up halfway between two samples and the parabolic peak is at 0.901229749985, which closely matches my answer. Parabolas are very good approximators for sinusoidal peaks, as in:

$$ \cos( x ) = 1 - \frac{x^2}{2} ... $$

Here is the code:

import numpy as np

#================================================
def main():

#---- Set Parameters

        p = 40
        n = 10

        phi = 0.0

#---- Construct the Averaged Signal

        x = np.zeros( 20 )

        omega = 2.0 * np.pi / p

        max_m = -1
        max_x = 0.0


        for m in range( 20 ):
          s = 0
          for k in range( n ):
            s += 1.0 * np.cos( omega * ( m - k ) + phi )

          x[m] = s / n

          print m, x[m]

          if max_x < x[m]:
             max_x = x[m]
             max_m = m

        print "Max:", max_m, max_x             

#---- Parabolic Approximation at Maximum

        vn = x[max_m-1]
        vz = x[max_m]
        vp = x[max_m+1]

        a = ( vn - 2.0 * vz + vp ) / 2
        b = ( vp - vn ) * 0.5
        c = vz

        d = -b / ( 2.0 * a )

        peak_m = max_m + d

        peak_x = a * d * d + b * d + c

        print d, peak_m, peak_x


# v = a x^2 + b x + c
# vn = a - b + c
# vz =         c
# vp = a + b + c

# a = ( vn - 2vz + vp ) / 2
# b = ( vp - vn ) / 2
# c = vz

# d = -b / (2a)          

#================================================
main()

---

Followup for MisterH:

This is the explanation for "calculation of bin 1 of a n-sized rectangular function within a p-sized DFT"

Here are your code lines:

DegreesPerSample<-360/p
sumxvalues<-sum(cos(((90+(0:(n-1))*DegreesPerSample))*pi/180)/n)
sumyvalues<-sum(sin(((90+(0:(n-1))*DegreesPerSample))*pi/180)/n)

Translated into math:

$$ dps = \frac{360}{p} $$

$$ \begin{aligned} S_x &= \sum_{s=0}^{n-1} \cos((90+s\cdot dps)\pi/180)/n \\ S_y &= \sum_{s=0}^{n-1} \sin((90+s\cdot dps)\pi/180)/n \\ \end{aligned} $$

Converting to radians.

$$ \beta = dps \cdot \frac{\pi}{180} = \frac{360}{p} \cdot \frac{\pi}{180} = \frac{2\pi}{p} $$

The units of $\beta$ are radians per sample.

Simplify $S_x$ using angle addition:

$$ \begin{aligned} S_x &= \frac{1}{n} \sum_{s=0}^{n-1} \cos(\pi/2 + \beta s ) \\ &= \frac{1}{n} \sum_{s=0}^{n-1} \left[\cos(\pi/2 )\cos(\beta s ) - \sin(\pi/2 )\sin(\beta s ) \right] \\ &= \frac{1}{n} \sum_{s=0}^{n-1} -\sin(\beta s ) \\ \end{aligned} $$

Similar for $S_y$.

$$ S_y \frac{1}{n} \sum_{s=0}^{n-1} \cos(\beta s ) $$

Introduce a "dummy" variable:

$$ k = 1 $$

$$ \begin{aligned} S_x[k] &= \frac{1}{n} \sum_{s=0}^{n-1} -\sin(\beta s k ) \\ S_y[k] &= \frac{1}{n} \sum_{s=0}^{n-1} \cos(\beta s k ) \\ \end{aligned} $$

Combine into a single complex equation:

$$ \begin{aligned} S[k] &= S_y[k] + i S_x[k] \\ &= \frac{1}{n} \sum_{s=0}^{n-1} \left[ \cos(\beta s k ) -i\sin(\beta s k ) \right] \\ &= \frac{1}{n} \sum_{s=0}^{n-1} e^{ -i \beta s k } \\ \end{aligned} $$

Define $ R_{0,n-1}[s] $ to be a unit rectangle function on the interval [0,n).

$$ \begin{aligned} S[k] &= \frac{1}{n} \sum_{s=0}^{n-1} e^{ -i \frac{2\pi}{p} s k } \\ &= \frac{1}{n}\left[ \sum_{s=0}^{p-1} R_{0,n-1}[s] \cdot e^{ -i \frac{2\pi}{p} s k } \right] \\ \end{aligned} $$

The variable names are a little different because I stuck to yours, but the expression in the brackets is the definition of the DFT applied to a rectangle function. When $k=1$ it matches your code.

Answer 2

Below is the analytic result for both the actual max value of $0.901243$ and the maximum value found by the OP of $0.898464$

The reason you are not getting the predicted maximum is your samples of the sine wave are not located exactly at the peak. This is clear if you zoom in on the plot and compare the two peak locations for the number of samples given (as I show in the plot below).

Also a direct approach to establish the frequency response for the moving average filter that properly accounts for both the phase and magnitude of the filter is as follows:

Consider the implementation of a moving average filter given as ($1/N$ scaling not shown):

The general transfer function for the moving average filter directly from the implementation block diagram (where $z^{-1}$ is the transfer function of a unit delay) is:

$$H(z) = \frac{1}{N}\sum_{n=0}^{N-1}z^{-n}$$

Using the well known relationship for the geometric series (see at bottom of answer) this is:

$$H(z) = \frac{1}{N}\frac{1-z^{-N}}{1-z}$$

And the frequency response (as the DTFT of the impulse response) is a continuous function of frequency $\omega$ for all $z = e^{j\omega}$, therefore the values of z on the unit circle. (Yes even though it is a discrete system the frequency response is indeed a continuous function and unique for $\omega$ over the range of $0$ to $\pi$ for real functions):

$$H(\omega) = \frac{1}{10}\frac{1-e^{-j\omega N}}{1-e^{-j\omega}}$$

$$=\frac{1}{N}\frac{e^{-j\omega N/2}(e^{+j\omega N/2} - e^{-j\omega N/2})}{e^{-j\omega /2}(e^{+j\omega /2} - e^{-j\omega /2})}$$

Using Euler's identity relating $\sin$ to the positive and negative exponential terms (see at bottom of answer) results in:

$$ H(\omega)= \frac{e^{j((N-1)/2)}}{N}\frac{sin(\omega N /2)}{sin(\omega/2)}$$

The exponential term has a magnitude of 1 for all $\omega$ but provides for the exact phase shift between input and output as evidenced in the plots.

The frequency $\omega$ is the normalized radian frequency in units of radians/sample, so in the OP's example $N=10$ and $\omega = \frac{2\pi}{40} = \frac{\pi}{20}$, and for$H(\pi/20)$ the result is:

$$H(\pi/20) = \frac{e^{j(4.5\pi/20)}}{10}\frac{sin(\pi/4)}{sin(\pi/40)} \approx 0.90124e^{-j0.70686}$$

Thus has a magnitude of $0.90124$ and an angle of $-.70686$ radians or -$40.500°$.

This result matches that provided by the freqz command in MATLAB/Octave:

>> h = freqz(ones(10,1),10, [0 pi/20]);
h = 1.0000 + 0.0000i  0.68351 - 0.58531i
>> abs(h(2))
ans = 0.90124
>> angle(h(2))
ans = -0.70686

Knowing the phase and amplitude from above, we can predict the result the OP got for the specific max sample point (15th sample) from $Asin(\omega n + \phi)$ using the amplitude $A$ and phase $\phi$ from above:

$$0.90124\sin(15 \pi/20 - 0.70686) = 0.898464$$

And the overall frequency response with magnitude in dB for $\omega = 0$ to $\pi$ is:

Note even though the samples of the output do not land on the exact peak does not mean this is not the peak of the output waveform. Consider Nyquists sampling theorem and how a waveform can be completely described (and if a single tone that would be its amplitude, frequency, phase etc) with relatively very few samples.

---

Relationships Used Above

Geometric Series

$$\sum_{n=0}^{N-1}r^k = \frac{1-r^N}{1-r}$$

Euler's Identity for sine

$$sin(\theta) = \frac{e^{+j\theta}-e^{-j\theta}}{2j}$$

Answer 3

The amplitude reduction is simply given as the magnitude of the transfer function of moving average filter. A moving average filter has a rectangular impulse response so the transfer function will be a $sinc()$ function. You need to sample the $sinc()$ function at the frequency or your sign wave

Answer 4

Puzzle solved, thanks to Cedron Dawg and Dan Boschen!

First, I ran a simple N point moving average of a sinewave, using the simulation model below:

I used the OP's values: N = 10, P = 40, sinewave amplitude = 1 and a simulation step size, $\Delta t$, equal to unity. The results, shown in the next figure, are the same as those of the OP:

The maximum amplitude of the (red) filtered sinewave was 0.8984644, same as the OP obtained. So far, so good.

Next, I ran the simulations with N = 100 for $\Delta t = 0.1$, N = 1000 for $\Delta t = 0.01$, and N = 10000 for $\Delta t = 0.001$. These gave maximum (red) filtered sinewave amplitudes of 0.9002978, 0.9003161, and 0.9003163, respectively. So sequentially reducing the step size by factors of ten, while increasing N by the corresponding factors of ten, results in the maximum filtered sinewave amplitude converging to approximately 0.9003163.

Now, if this was an analog system, with continuous time, the magnitude of the transfer function would be

$$|H(\omega)| = \frac{\tau_a}{\tau_i} \times \operatorname{sinc}(f\tau_a) = \frac{\tau_a}{\tau_i} \times \frac{\sin(\pi f\tau_a)}{\pi f\tau_a} \tag{1}$$

where $\tau_a = 10$ is the integration aperture (aka 'gate'), $\tau_i = 10$ is the integration time constant, and f = 1/P = 1/40 is the frequency. As per Hilmer, the impulse response is simply a rectangular pulse with amplitude = $1/\tau_i$, duration = $\tau_a$ and $\tau_a = \tau_i$ for unity gain as an averager. With $\tau_a$, $\tau_i$ and f substituted into equation (1), the result is

$$|H(\omega)| = 0.9003163162 \tag{2}$$

Answer 5

From a slightly less "dsp-like" point of view, slightly more geometric / time series, but this also works:

The relation between the sinusoid (of amplitude 1) and the unit circle is well known.

Instead of thinking of a moving average as a geometric mean on a window that slides from left to right over the time series, you could also define it as the cumulative sum of a lower ($*1/n$) amplitude sine (from right to left in the window), as the window slides from left to right over the time series: I refer to the left panel in the image below:

(best to open the image in a new tab)

Now looking at the right panel in the image above, at the unit circle: the average can be seen as the cumulative sum of the vectors formed between the circle center and the purple dots on the small circle: the x- and y-values are respectively the cos & sin of $90° + 0:(n-1) * 360/period * 1/n$). This cumulative vector sum also lies on a circle: the dotted cyan circle's centre lies at $(0.6353,0.05002)$, and has a radius of $0.63726$. The horizontal coordinate of the circle center lies at +/- $1/(2*n)$.

You then calculate the length from the end point of the circular segment to the origin, using the sum of the x-values and y-values via pythagoras' theorem, and you get your amplitude reduction, as indicated by the length of the radius of the dotted red circle, in this case $0.9012426$. The circle will rotate if you chose a different start point (not $90°$), but the end of the circular segment will always be on the same dotted red circle.

DegreesPerSample<-360/p
sumxvalues<-sum(cos(((90+(0:(n-1))*DegreesPerSample))*pi/180)/n)
sumyvalues<-sum(sin(((90+(0:(n-1))*DegreesPerSample))*pi/180)/n)
((sumxvalues^2)+(sumyvalues^2))^.5

I believe there are parallels with the Hilbert Transform's In-Phase and Quadrature components?

The lag of a simple moving average is $(n-1)/2$ (I believe you dsp guys call this group delay). This is also the center of gravity of the rectangular impulse response. You can also read that off the circular phase delay plot: for $p=40$ and $n=10$, the cyan arrow shows that the angle of the end of the circular segment to the origin is $130.5$ degrees: $90$ (the startpoint)$ + (n-1)/2 * 360/p$. It can also be seen as $4.5$: $(n-1)/2$ small purple circle segments on the unit circle.

But there is more information to be found here:

there are causal filters with negative weights at the back (left side) of the window. Using negative weights, you can create a causal filter that is "in-phase" with a sinusoidal signal. In the $n=10, p=40$ case, instead of the (SMA) weights:

0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

if you use:

-0.1000 -0.1000 -0.1000  0.0764  0.1000  0.1000  0.1000  0.1000  0.1000  0.1000

which is the equivalent of $31.17959$ % negative weights, the adapted simple moving average filter's result will be "in-phase" with the sinusoid.

If you use the same amplitude reduction concept, and you multiply this new filter's output's amplitude with $2.051051$, you get your original sinusoid, constructed from a causal filter on a window of length $n$. A quick calculation learns that the sum of the weights of the new amplitude-adjusted weights is $0.7720322$.

All good and well when you know the period of your sinusoid. But how could you create the so-called xvalues (cosines of the angles) from the yvalues (values of "a" signal in the window $/n$)? For an arbitrary -unknown- period and amplitude?

---

Updated 10/04/20:

So, how do you calculate the % of negative weights at the left side of the window in the convolution such that the filter's result is in-phase with a sinusoid of period p? You have 2 percentages: the negative (p1), and the other one, so I called it "p1function":

p1function<-function(n,DegreesPerSample)
{
  xvalues3<-round(cos((90+((0:(n-1))*DegreesPerSample))*pi/180)*smavector(n),5)
  yvalues3<-round(sin((90+((0:(n-1))*DegreesPerSample))*pi/180)*smavector(n),5)
  outputc<-circleFromThreePoints(cumsum(xvalues3)[1],cumsum(xvalues3)[2],cumsum(xvalues3)[3],cumsum(yvalues3)[1],cumsum(yvalues3)[2],cumsum(yvalues3)[3])
  hc<-outputc[1];kc<-outputc[2];rc<-outputc[3];
  xc<-seq((hc-rc),(hc+rc),length.out=1001)
  x1<-hc
  y1<-kc
  x2<-(sum(xvalues3)/2)
  y2<-kc+(((rc^2)-(((sum(xvalues3)/2)-hc)^2))^0.5)
  x3<-0
  y3<-sinn(90+DegreesPerSample)/n
  x4<-sum(xvalues3)
  y4<-sum(yvalues3)
  startangle<-atan2((y3-y1),(x3-x1))*180/pi
  midangle<-atan2((y2-y1),(x2-x1))*180/pi
  endangle<-atan2((y4-y1),(x4-x1))*180/pi
  endangle<-ifelse(endangle<0,endangle+360,endangle)
  p1est<-100*((endangle+startangle)-(midangle+startangle))/(endangle+startangle)
  return(p1est)
}

Which uses another function to estimate a circle from 3 points:

circleFromThreePoints<-function(x1,x2,x3,y1,y2,y3)
{
  vara<-x1*(y2-y3)-y1*(x2-x3)+x2*y3-x3*y2;
  varb<-(x1*x1+y1*y1)*(y3-y2)+(x2*x2+y2*y2)*(y1-y3)+(x3*x3+y3*y3)*(y2-y1);
  varc<-(x1*x1+y1*y1)*(x2-x3)+(x2*x2+y2*y2)*(x3-x1)+(x3*x3+y3*y3)*(x1-x2);
  vard<-(x1*x1+y1*y1)*(x3*y2-x2*y3)+(x2*x2+y2*y2)*(x1*y3-x3*y1)+(x3*x3+y3*y3)*(x2*y1-x1*y2)
  varx<- -varb/(2*vara)
  vary<- -varc/(2*vara)
  varr<- (((varb*varb)+(varc*varc)-(4*vara*vard))/(4*vara*vara))^0.5
  # x, y , r: 
  # (x-x1)^2+(y-y1)^2 = r^2 
  # h,k,r for equation: (x-h)^2+(y-k)^2 = r^2
  # To plot: upp<-(((r^2)-((x-h)^2))^0.5)+k & dwn<--(((r^2)-((x-h)^2))^0.5)+k
  return(c(round(varx,5),round(vary,5),round(varr,5)))
}

The hashtag means it's a comment.

And of course this filter's output's amplitude is not the same, so, going back to the plot I made above, you just need to adjust its height: same way as above:

ampfactor<-function(n,p1est,DegreesPerSample)
{
  xvalues3<-cos((90+((0:(n-1))*DegreesPerSample))*pi/180)/n
  yvalues3<-sin((90+((0:(n-1))*DegreesPerSample))*pi/180)/n
  1/(sum(xvalues3*onesfunc(n,p1est))^2+sum(yvalues3*onesfunc(n,p1est))^2)^0.5
}

The fun thing is, I think somewhere in here there's a causal version of my favorite Hodrick-Prescott Filter. Could lead to an accurate instantaneous frequency estimator, as it only needs 3 points -given the circle-, no?

As one varies the p1%, you get closer to the wave. This works well in theory. I do believe using negative weights at the back of the window can bring you just 1 more step closer. Last step would be using this concept on random data. Obviously causality cannot be broken, but you can get closer. Please do correct me where needed.