# Next steps for analyzing this PWM algorithm and choosing a low pass filtering technique?

The Broadcom BCM2711 GPIO controller used by the Raspberry Pi has a hardware PWM implementation option using an algorithm that provides an average value defined by a rational fraction $$N/M$$ defined as follows in section 8.3. PWM Implementation starting on page 156:

1. Set context = 0
2. context = context + N
3. if (context >= M) context = context - M send 1 else send 0
4. Repeat from step 2

where context is a register which stores the result of the additions/subtractions.

I've implemented this in a Python script to try to start to get an idea how I'll need to filter it in order to get a low-noise, stable voltage. Here with $$M=200$$ I've incremented through $$0 \le N \le M$$ for over 1000 time points each, then applied a gaussian smoothing with a $$\sigma$$ of 10 time points just to get a quick look.

I took the Fourier transform of the middle 1024 points, ignoring a $$3 \sigma$$ buffer on each side.

I was surprised to see that some values of $$N$$ produce signals that are orders of magnitude more noisy than other nearby values of $$N$$.

Questions:

1. This kind of analysis is new territory for me, if I've made any major errors in the analysis I would appreciate knowing about it.
2. If I need near-DC voltages I can just run the BCM2711 at a moderately high frequency and filter the heck out of the output with a simple RC filter. But if I wanted to work closer to the frequency that the algorithm was running, what additional analysis should I do on this algorithm before deciding on a more aggressive low-pass filter and simulating it's performance?

Here's the Python script for the plots and analysis:

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
from numpy.fft import rfft, rfftfreq

sigma = 10
ignore = 3 * sigma
M = 200
big = []
for N in range(M+1):
context = 0
bob = []
for i in range(1024 + 2 * ignore):
context += N
if (context >= M):
context -= M
bob.append(1)
else:
bob.append(0)
big.append(bob)

big = np.array(big, dtype=float)
fbig = ndimage.gaussian_filter1d(big, sigma=sigma, mode='mirror', order=0, axis=-1)
fbig_trunc = fbig[:, ignore:-ignore]
print('fbig_trunc.shape: ', fbig_trunc.shape)
ft = rfft(fbig_trunc, axis=-1)
power = np.abs(ft)**2
power /= power.max()
dB_rel = 10 * np.log10(power)
print('power.shape: ', power.shape)
print('power.max(): ', power.max())
print('dB_rel.shape: ', dB_rel.shape)
print('dB_rel.max(): ', dB_rel.max())

if True:
plt.figure()
for thing in fbig:
plt.plot(thing, linewidth=0.8)
plt.ylim(-0.02, 1.02)
plt.title('sigma={:2d}, M={:3d}'.format(sigma, M), fontsize=16)
plt.xlabel('time steps', fontsize=16)
plt.show()

if True:
plt.figure()
maxi = dB_rel[:, 1:201].max()
plt.imshow(dB_rel[:, :201], vmax=maxi, vmin=-100)
plt.colorbar()
plt.gca().set_aspect(2.0)
plt.title('FT relative power (dB), sigma={:2d}, M={:3d}'.format(sigma, M), fontsize=16)
plt.xlabel('frequency', fontsize=16)
plt.ylabel('N', fontsize=16)
plt.show()

if True:
plt.figure()
std = fbig_trunc.std(axis=-1)
plt.plot(10*np.log10(std**2))
plt.title('stdev relative power, sigma={:2d}, M={:3d}'.format(sigma, M), fontsize=16)
plt.xlabel('N', fontsize=16)
plt.ylabel('stdev (dB rel)', fontsize=16)
plt.show()