2
\$\begingroup\$

I've written (and tested) a simple least mean square adaptive filter . The steps of the algorithm are:

$$ For n = 0, 1, 2, \ldots: $$
    $$ x = [u(n), u(n+1), \ldots, n(n-p+1)]^T $$
    $$ e(n) = d(n) - w(n)^Tx(n) $$
    $$ w(n+1) = w(n) + \mu*e(n)^Tx(n) $$

Where \$u\$ is the input signal, \$w\$ are the weights of the filter, \$p\$ is the order (number of taps), \$e\$ is the error signal, and \$d\$ is the desired signal.

import numpy as np
import matplotlib.pyplot as plt
import traceback
import sys

class LMS:
    def __init__(self, input_signal, desired_signal, num_taps, learning_rate):
        self.u = input_signal
        self.d = desired_signal
        self.num_taps = num_taps
        self.mu = learning_rate
        self.num_points = len(self.u)
        self.weights = np.zeros(num_taps)
        self.y = np.zeros(self.num_points)
        self.e = np.zeros(self.num_points)

    def equalize(self):
        for n in xrange(self.num_taps, self.num_points):
            x = self.u[n:n-num_taps:-1]
            self.y[n] = np.dot(x, self.weights)
            self.e[n]= self.d[n] - self.y[n]
            self.weights = self.weights + self.mu*x*self.e[n] 

try:
    np.random.seed(1337)
    ulen    = 20000
    coeff   = np.concatenate(([1], np.zeros(10), [-0.9], np.zeros(7), [0.1]))
    u       = np.random.randn(ulen)
    d       = np.convolve(u, coeff)
    num_taps = 20 
    step = 0.01 

    eq = LMS(u, d, num_taps, step)
    eq.equalize()
    plt.figure()
    plt.semilogy()
    plt.subplot(1,1,1)
    plt.plot(np.abs(eq.e))
    plt.show()
except:
    type, value, tb = sys.exc_info()
    traceback.print_exc()
    pdb.post_mortem(tb)

As written the algorithm converges (I get the following plot of the error vector):

Error vs iteration number

\$\endgroup\$
3
\$\begingroup\$

Exceptions

You don't need the sys.exc_info call since pdb.post_mortem without argument uses the exception being currently handled. You don't need pdb either since it serves for development purposes and, at this stage, you should have working code already.

And since all what's left in the except clause is traceback.print_exc(), you could remove the try .. except and let the interpreter print the traceback for you.

Also, it's a very bad habbit to use bare excepts. You should know what kind of exceptions you’re willing to handle. (Even though it can be usefull with pdb to examine various causes of issues.)

Class

Given how you use the class (initialization of attributes, computation and getting one attribute back), I think you could do equally well with a simple function. Just initialize what you need before the for loop and return the computed array:

def LMS(input_signal, desired_signal, num_taps, learning_rate):
    num_points = len(input_signal)
    weights = np.zeros(num_taps)
    equalized = np.zeros(num_points)

    for n in xrange(num_taps, num_points):
        x = input_signal[n:n-num_taps:-1]
        equalized[n]= desired_signal[n] - np.dot(x, weights)
        weights = weights + learning_rate*x*equalized[n]

    return equalized

Note that I kept the meaningfull names to help read the algorithm. I also removed y which was a temporary variable and served no purpose after computation.

main

It is also good practice to avoid keeping code at the top-level of the file. Just in case you want to import your file (for testing purposes in an interactive shell, for instance), you should guard the code in an if __name__ == '__main__': clause.

The whole code becomes:

import numpy as np
import matplotlib.pyplot as plt

def LMS(input_signal, desired_signal, num_taps, learning_rate):
    num_points = len(input_signal)
    weights = np.zeros(num_taps)
    equalized = np.zeros(num_points)

    for n in xrange(num_taps, num_points):
        x = input_signal[n:n-num_taps:-1]
        equalized[n]= desired_signal[n] - np.dot(x, weights)
        weights = weights + learning_rate*x*equalized[n]

    return equalized


if __name__ == '__main__':
    np.random.seed(1337)
    ulen = 20000
    coeff = np.concatenate(([1], np.zeros(10), [-0.9], np.zeros(7), [0.1]))
    u = np.random.randn(ulen)
    d = np.convolve(u, coeff)
    num_taps = 20 
    step = 0.01 

    eq = LMS(u, d, num_taps, step)

    plt.figure()
    plt.semilogy()
    plt.subplot(1,1,1)
    plt.plot(np.abs(eq))
    plt.show()

(I removed the extra whitespace arround assignements since it is recommended by PEP8)

\$\endgroup\$
  • \$\begingroup\$ Thank you! Ironically, I refactored this from some code that was just a function. I want it as a class, as it's part of a project with several types of equaliser (which will decide the eq type at runtime). \$\endgroup\$ – Tom Kealy Dec 8 '15 at 0:14
  • 1
    \$\begingroup\$ @TomKealy You should have mentionned that in the question ;) Anyway, for such project, you probably want to invert your logic: call your class Equalizer and each method beside the constructor will be named after the type of operation it performs. \$\endgroup\$ – 409_Conflict Dec 8 '15 at 9:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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