HiveBrain v1.2.0
Get Started
← Back to all entries
patternpythonMinor

Least Mean Square channel equalizer

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
meanequalizerleastsquarechannel

Problem

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):

Solution

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)

Code Snippets

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
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()

Context

StackExchange Code Review Q#113164, answer score: 3

Revisions (0)

No revisions yet.