# Recursive Least Squares (RLS)¶

New in version 0.1.

Changed in version 1.0.0.

The Recursive Least Squares filter [1] can be created as follows

>>> import padasip as pa
>>> pa.filters.FilterRLS(n)


where the n is amount of filter inputs (size of input vector).

## Algorithm Explanation¶

The RLS adaptive filter may be described as

$$y(k) = w_1 \cdot x_{1}(k) + ... + w_n \cdot x_{n}(k)$$,

or in a vector form

$$y(k) = \textbf{x}^T(k) \textbf{w}(k)$$,

where $$k$$ is discrete time index, $$(.)^T$$ denotes the transposition, $$y(k)$$ is filtered signal, $$\textbf{w}$$ is vector of filter adaptive parameters and $$\textbf{x}$$ is input vector (for a filter of size $$n$$) as follows

$$\textbf{x}(k) = [x_1(k), ..., x_n(k)]$$.

The update is done as folows

$$\textbf{w}(k+1) = \textbf{w}(k) + \Delta \textbf{w}(k)$$

where $$\Delta \textbf{w}(k)$$ is obtained as follows

$$\Delta \textbf{w}(k) = \textbf{R}(k) \textbf{x}(k) e(k)$$,

where $$e(k)$$ is error and it is estimated according to filter output and desired value $$d(k)$$ as follows

$$e(k) = d(k) - y(k)$$

The $$\textbf{R}(k)$$ is inverse of autocorrelation matrix and it is calculated as follows

$$\textbf{R}(k) = \frac{1}{\mu}( \textbf{R}(k-1) - \frac{\textbf{R}(k-1)\textbf{x}(k) \textbf{x}(k)^{T} \textbf{R}(k-1)} {\mu + \textbf{x}(k)^{T}\textbf{R}(k-1)\textbf{x}(k)} )$$.

The initial value of autocorrelation matrix should be set to

$$\textbf{R}(0) = \frac{1}{\delta} \textbf{I}$$,

where $$\textbf{I}$$ is identity matrix and $$\delta$$ is small positive constant.

## Stability and Optimal Performance¶

Make the RLS working correctly with a real data can be tricky. The forgetting factor $$\mu$$ should be in range from 0 to 1. But in a lot of cases it works only with values close to 1 (for example something like 0.99).

## Minimal Working Examples¶

If you have measured data you may filter it as follows

import numpy as np
import matplotlib.pylab as plt

# creation of data
N = 500
x = np.random.normal(0, 1, (N, 4)) # input matrix
v = np.random.normal(0, 0.1, N) # noise
d = 2*x[:,0] + 0.1*x[:,1] - 4*x[:,2] + 0.5*x[:,3] + v # target

# identification
f = pa.filters.FilterRLS(n=4, mu=0.1, w="random")
y, e, w = f.run(d, x)

# show results
plt.figure(figsize=(15,9))
plt.plot(d,"b", label="d - target")
plt.plot(y,"g", label="y - output");plt.legend()
plt.subplot(212);plt.title("Filter error");plt.xlabel("samples - k")
plt.plot(10*np.log10(e**2),"r", label="e - error [dB]");plt.legend()
plt.tight_layout()
plt.show()


An example how to filter data measured in real-time

import numpy as np
import matplotlib.pylab as plt

# these two function supplement your online measurment
def measure_x():
# it produces input vector of size 3
x = np.random.random(3)
return x

def measure_d(x):
# meausure system output
d = 2*x[0] + 1*x[1] - 1.5*x[2]
return d

N = 100
log_d = np.zeros(N)
log_y = np.zeros(N)
filt = pa.filters.FilterRLS(3, mu=0.5)
for k in range(N):
# measure input
x = measure_x()
# predict new value
y = filt.predict(x)
# do the important stuff with prediction output
pass
# measure output
d = measure_d(x)
# update filter
# log values
log_d[k] = d
log_y[k] = y

### show results
plt.figure(figsize=(15,9))
plt.plot(log_d,"b", label="d - target")
plt.plot(log_y,"g", label="y - output");plt.legend()
plt.subplot(212);plt.title("Filter error");plt.xlabel("samples - k")
plt.plot(10*np.log10((log_d-log_y)**2),"r", label="e - error [dB]")
plt.legend(); plt.tight_layout(); plt.show()


## References¶

 [1] Ali H Sayed and Thomas Kailath. Recursive least-squares adaptive filters. The Digital Signal Processing Handbook, pages 21–1, 1998.

## Code Explanation¶

class padasip.filters.rls.FilterRLS(n, mu=0.99, eps=0.1, w='random')[source]

Bases: padasip.filters.base_filter.AdaptiveFilter

Args:

• n : length of filter (integer) - how many input is input array (row of input matrix)

Kwargs:

• mu : forgetting factor (float). It is introduced to give exponentially less weight to older error samples. It is usually chosen between 0.98 and 1.

• eps : initialisation value (float). It is usually chosen between 0.1 and 1.

• w : initial weights of filter. Possible values are:

• array with initial weights (1 dimensional array) of filter size
• “random” : create random weights
• “zeros” : create zero value weights
adapt(d, x)[source]

Adapt weights according one desired value and its input.

Args:

• d : desired value (float)
• x : input array (1-dimensional array)
run(d, x)[source]

This function filters multiple samples in a row.

Args:

• d : desired value (1 dimensional array)

• x
: input matrix (2-dimensional array). Rows are samples,

columns are input arrays.

Returns:

• y : output value (1 dimensional array). The size corresponds with the desired value.
• e : filter error for every sample (1 dimensional array). The size corresponds with the desired value.
• w : history of all weights (2 dimensional array). Every row is set of the weights for given sample.