3. Signal Statistics and Noise

3.1 Moving Average

We can use a circular buffer (an array and a moving index) to store the last N samples of a signal. The moving average $\text{MA}$ using a window size of $N$ can be expressed as:

$\text{MA}(t) = \frac{1}{N} \sum_{i=0}^{N-1} x(t-i)$

The sum holds $\sum_{i=0}^{N-1} x(t-i)$ , and average stores $\text{MA}(t)$ . The value is updated efficiently as follows:

$\text{New Sum} = \text{Old Sum} - x(t - N) + x(t)$

This allows you to then update $\text{MA}(t)$ as:

$\text{MA}(t) = \frac{\text{New Sum}}{N}$

In code, this looks like:

static float32_t circular_buffer[WINDOW_SIZE];
static int circular_index = 0;  // Index where the next element will go into the circular buffer
static float32_t sum = 0.0;     // Keep track of the sum to calculate the moving average

static void print_signal_and_moving_average_efficient(void)
{
    for (int i = 0; i < KHZ1_15_SIG_LEN; i++)
    {
        // Remove the oldest element from the sum
        sum -= circular_buffer[circular_index];

        // Add the new element to the circular buffer and to the sum
        circular_buffer[circular_index] = input_signal_f32_1kHz_15kHz[i];
        sum += circular_buffer[circular_index];

        // Compute and print the moving average if the window is filled
        if (i >= WINDOW_SIZE - 1)
        {
            float32_t average = sum / WINDOW_SIZE;
            printf("%f, %f\r\n", input_signal_f32_1kHz_15kHz[i], average);
        }

        // Move the index, wrapping around if necessary
        circular_index = (circular_index + 1) % WINDOW_SIZE;
    }
}

The result can be visualized using the Serial Plotter:

mean_result

Where the blue line is the input signal, and the red line is the moving average. Also see that there is a delay of $N$ samples between the input signal and the moving average. There are also discontinuities at the beginning and end of the original signal.

3.2 Moving Standard Deviation

We can calculate the moving standard deviation efficiently by keeping track of both the sum of the elements and the sum of their squares in the window. Here's the formula to calculate the standard deviation for a window of $N$ elements:

$\text{StdDev}(t) = \sqrt{\frac{1}{N} \sum_{i=0}^{N-1} x(t-i)^2 - \left(\frac{1}{N} \sum_{i=0}^{N-1} x(t-i)\right)^2}$

This can be simplified to:

$\text{StdDev}(t) = \sqrt{\frac{1}{N} \text{sumOfSquares} - \left(\frac{1}{N} \text{sum}\right)^2}$

In code, this looks like:

#include <math.h> // for sqrt()

#define WINDOW_SIZE 20

static float32_t circular_buffer[WINDOW_SIZE];
static float32_t sum_of_squares = 0.0;
static int circular_index = 0;
static float32_t sum = 0.0;

static void print_signal_and_moving_std_dev(void)
{
    for (int i = 0; i < KHZ1_15_SIG_LEN; i++)
    {
        float32_t old_value = circular_buffer[circular_index];
        float32_t new_value = input_signal_f32_1kHz_15kHz[i];

        // Update the sum and sum of squares
        sum -= old_value;
        sum += new_value;

        sum_of_squares -= old_value * old_value;
        sum_of_squares += new_value * new_value;

        // Replace the oldest value in the circular buffer
        circular_buffer[circular_index] = new_value;

        if (i >= WINDOW_SIZE - 1)
        {
            float32_t mean = sum / WINDOW_SIZE;
            float32_t std_dev = sqrt((sum_of_squares / WINDOW_SIZE) - (mean * mean));
            printf("%f, %f\r\n", new_value, std_dev);
        }

        // Move the index
        circular_index = (circular_index + 1) % WINDOW_SIZE;
    }
}

The result can be visualized using the Serial Plotter:

std_dev_result

3.3 Speeding Up the Moving Standard Deviation

We can speed up the moving standard deviation by using the arm_sqrt_f32 function from the CMSIS-DSP library. It is designed for high performance on ARM processors, especially those with a Floating-Point Unit (FPU).