Implementing logarithmic AGC (automatic gain control)

Question

I'd like to implement the Log(arithmic) AGC from

• Understanding Digital Signal Processing

But with

import numpy as np
import matplotlib.pyplot as plt

def signal():
    cycles = 15
    resolution = cycles * 20
    length = 2 * np.pi * cycles
    x = np.cos(np.arange(0, length, length / resolution))
    cycles = 10
    resolution = cycles * 20
    length = 2 * np.pi * cycles
    x = np.append(x, [2.5 * np.cos(np.arange(0, length, length / resolution))])
    cycles = 15
    resolution = cycles * 20
    length = 2 * np.pi * cycles
    x = np.append(x, [np.cos(np.arange(0, length, length / resolution))])
    cycles = 20
    resolution = cycles * 20
    length = 2 * np.pi * cycles
    x = np.append(x, [0.25 * np.cos(np.arange(0, length, length / resolution))])
    return x

def logarithmic_agc_book(x, alpha, R):
    g1 = 0
    y = np.zeros((len(x),))
    for n, xn in enumerate(x):
        y[n] = xn * np.exp2(g1)
        tmp = y[n] * y[n]
        # tmp = LPF(tmp)
        tmp = 2 * np.log2(R) - np.log2(tmp)
        tmp = tmp * alpha
        g = tmp + g1
        g1 = g
    return y


def logarithmic_agc_web(x, alpha, R):
    g1 = 0
    y = np.zeros((len(x),))
    for n, xn in enumerate(x):
        y[n] = xn * np.exp(g1)
        tmp = np.abs(y[n])
        tmp = np.log(R) - np.log(tmp)
        tmp = tmp * alpha
        g = tmp + g1
        g1 = g
    return y


x = signal()
fig1, (ax1, ax2) = plt.subplots(nrows=2, ncols=1)
fig1.tight_layout()
ax1.plot(x, label="Signal")
ax1.set_xlim([0, 1199])
ax1.set_ylim([-5, 5])
ax1.set_xlabel(f"Time (samples)")
ax1.legend()
ax2.plot(logarithmic_agc_book(x, 0.01, 1), label="Logarithmic AGC (book)")
ax2.plot(logarithmic_agc_web(x, 0.01, 1), label="Logarithmic AGC (web)")
ax2.set_xlim([0, 1199])
# ax2.set_ylim([-5, 5])
ax2.set_xlabel(f"Time (samples)")
ax2.legend()
plt.show()

I only get for both variants

which seems that both are not correct (e.g., peaks at ~50). Does somebody spot the problem in the implementations?

Answer 1

I believe the issue is in deriving the power with the absolute value or squaring only, which when the log of this is taken, results in very large negative numbers as the result approaches zero in each case. When negated from the target, this produces very large peak values as the error function. The loop will then average this to zero error. The large peaks create a significant offset in this overall average.

We see this in the plots below for the OP's AGC (web) test case showing the Error Signal (result of np.og(R)-np.log(tmp)), and log of the output (result of np.log(tmp).

Using either an absolute value or squaring function is ok to do, but a moving average or other simple low pass filter is required at the output prior to taking the log to complete the power detector operation appropriately. It is important that the low pass filter have a cutoff frequency that is significantly higher than the loop bandwidth to not interfere with the loop (so relatively short averaging time). The use of a squaring function results in a "true-rms" power detector, while with the absolute value we still get a metric proportional to output power but it is not "true-rms". We have the same option when choosing hardware performing the power detector function in analog and microwave systems.

A non-true rms power detector will result in an output power level that is dependent on the peak-avg of the waveform, so is therefore waveform dependent. When this is a concern, a true-rms power detector should be used.

The average of this error signal is zero, which is what would be expected to happen for this first order loop as implemented.

Zooming in provides insight into what is happening as we approach abs(0) prior to the log:

And this plot is showing the log of the absolute value of the output directly. The solution is to do a simple moving average or other low pass filter on the absolute value output prior to taking the log.

I updated the AGC (web) case with a simple 10 sample moving average using the python collections.deque function as a FIFO:

from collections import deque
def logarithmic_agc_web(x, alpha, R):
    g1 = 0
    avg_dur=10
    mavg = deque(np.zeros(avg_dur),avg_dur)
    y = np.zeros((len(x),))
    for n, xn in enumerate(x):
        y[n] = xn * np.exp(g1)
        mavg.append(np.abs(y[n]))
        tmp = np.sum(mavg)/avg_dur
        tmp = alpha * (np.log(R) - np.log(tmp))
        g1 = tmp + g1
    return y

Note that the average value of a sine wave is 0.637 times the peak value. The resulting output is in good agreement with expected $1/0.637 \approx 1.57$ as we see zooming in on the magnitude of the output below:

Please also see this for a useful and efficient computation of the logs for this AGC implementation (and why implementing with a log2 instead of a natural log is preferred): Efficient Log2 and dB from Floating Point and Fixed Point Representation

Note: This AGC (with expanded details on why a log and antilog function is necessary, and how to model the loop and compute the loop bandwidth, what bandwidth should be used, etc...), many other key functional building blocks for Software Radio, and tons of examples in Python are featured in an online course I have starting this week- more info here: "DSP For Software Radio"