6
$\begingroup$

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

  • Understanding Digital Signal Processing enter image description here

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

enter image description here

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

$\endgroup$
3
  • 2
    $\begingroup$ I can't help you with numpy. But I can tell you that the multiplier on the right side of my Figure 13-77 computes power, not magnitude. $\endgroup$ Commented Sep 19, 2023 at 9:25
  • $\begingroup$ Thanks, I will update the figure. $\endgroup$
    – ge45mue
    Commented Sep 19, 2023 at 9:28
  • $\begingroup$ BTW, you might wanna consider feedforward AGC, because you can put the envelope detection (and all of the AGC math) on the undelayed input signal and then apply the net gain change to a slightly delayed signal and start adjusting gain before the sudden amplitude change occurs. $\endgroup$ Commented Sep 19, 2023 at 18:08

1 Answer 1

6
$\begingroup$

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.

error signal

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:

enter image description here

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.

ln(abs(y))

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

correct result

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:

enter image description here

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 that runs a couple times each year - more info here on the latest course offerings: https://dsprelated.com/courses and https://ieeeboston.org/courses/

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.