# Amplitude reduction of Simple Moving Average on sinusoid

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".

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.

$$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.

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"

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}

$$\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.

• Excellent! Thanks for solving the puzzle! – Ed V Mar 31 at 13:41
• @EdV You're welcome. Since the peak occurs halfway between two samples, the original signal can be shifted by a half sample (phi = np.pi/p) and the program then finds the true peak value which matches the results of my original analysis. Frequency response is generally done with $e^{i\omega}$ (like Dan did) instead of the signal definition, but I wanted a more concrete demonstration. It's great when the signal definition reemerges at the end after that mess. – Cedron Dawg Mar 31 at 14:04
• That was a beautiful analysis you did and I upvoted already! One last thing: should I delete my "answer"? I don't want to crud the place up with deadwood. – Ed V Mar 31 at 14:07
• Thank you. Absolutely not! It is more illustrative than the original question. I am going to edit your Tex just a tad. – Cedron Dawg Mar 31 at 14:10
• @CedronDawg Yes you were right about the 0.90124. I actually came up with another way to solve it today (Covid-19 is good for creativity), which returns the same result (in R): sumxvalues<-sum(cos(((90+(0:(n-1))*(360/p)))*pi/180)/n) sumyvalues<-sum(sin(((90+(0:(n-1))*(360/p)))*pi/180)/n) ((sumxvalues^2)+(sumyvalues^2))^.5 which for p=40,n=10 returns 0.9012426 If you like I can elaborate how I got there? – MisterH Apr 2 at 20:43

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}$$

• I upvoted this earlier and now you have made it even more thorough! Beautiful work and I would give it another upvote if I could! – Ed V Mar 31 at 22:38
• Oh thank you! Yes much more concise now than what I first had. Always fun to play with math and it keeps us sharp, right? – Dan Boschen Mar 31 at 22:41
• @DanBoschen sorry about the 0.901243 thing: 0.898464 was an incorrect value for the maximum amplitude value of the filtered result. – MisterH Apr 2 at 20:49
• @MisterH no need to apologize! It is the maximum value at the sample location you had (the output simply wasn't sampled right at its max value). – Dan Boschen Apr 2 at 20:55

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

• This is incorrect. See my side note. – Cedron Dawg Mar 30 at 23:07

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}$$

• A puzzle indeed. Your analysis looks solid and I don't see a mistake in my answer either. I've done this math before in different form so the answer I got is what I expected. See the "As N gets large" section in dsprelated.com/showarticle/1038.php or my answer in dsp.stackexchange.com/questions/63076/… I'm taking a closer look at this. – Cedron Dawg Mar 31 at 3:18
• Discrepancy solved, see my followup. – Cedron Dawg Mar 31 at 3:58

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.

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),cumsum(xvalues3),cumsum(xvalues3),cumsum(yvalues3),cumsum(yvalues3),cumsum(yvalues3))
hc<-outputc;kc<-outputc;rc<-outputc;
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.

• I'm checking it out. This may take more time than I have right now, so please be patient for a reply. First impression, your variable name "phaseshift" should really be "DegreesPerSample", "dps", freq_dps, or similar, as it is a frequency parameter. I would gently suggest that you transition from degrees to radians. It will simplify your expressions some and make it more palatable for some. $$dps \left( \frac{Degrees}{Sample} \right) = \frac{ 360\left( \frac{Degrees}{Cycle} \right) }{ p \left( \frac{Samples}{Cycle} \right) }$$ – Cedron Dawg Apr 3 at 17:52
• NIcely done graphic, BTW. – Cedron Dawg Apr 3 at 17:58
• @CedronDawg no time constraint, I appreciate your input. And yes, radians undoubtedly makes more sense in a signal processing thread. – MisterH Apr 5 at 14:34
• Sorry I don't have time to spell it out, so I leave it as a challenge. Unfurled and such your code example can be framed as the equivalent of a calculation of bin 1 of a n-sized rectangular function within a p-sized DFT, thus the alias sinc function emerges. I'm still not getting how you got there, but it's cool and I'll get there eventually. – Cedron Dawg Apr 5 at 16:19
• @CedronDawg You're a generous person.It will take me some time on wolframalpha to understand your equations; there ought to be a relatively simple way to redefine my "xvalues3" variable for non-sinusoidal, quasi-random univariate data. It's a bit frustrating right now, but thank you very much! – MisterH Apr 11 at 22:02