Feeds:
Posts

## Matlab and Fibonacci Numbers

Matlab has long been an immensely convenient tool for processing and visualizing signals in a DSP system. The tool makes it very easy to simultaneously visualize a signal in multiple domains. The following is a simple demonstration of some of these signal-processing capabilities.

The Fibonacci sequence is the sequence of numbers formed when each number in the sequence is equal to the sum of the two Fibonacci numbers immediately preceding that number. The sequence is initialized by the numbers 0 and 1, and looks something like this: 0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,etc. It is trivial to design a digital filter that produces the Fibonacci sequence as it’s impulse response. The impulse response of a discrete-time or digital filter is the response produced when the filter is stimulated by the Kronecker delta function:

$\delta[n] = \begin{cases} 1, & n = 0 \\ 0, & otherwise \end{cases}$

If the initial state of such a filter is such that all internal elements are at zero, the filter described by the following difference equation can produce an impulse response equal to the Fibonacci sequence. The input and output sequences are denoted x[n] and y[n], respectively.

$y[n] = x[n] + y[n-1] + y[n-2]$

Taking the Z-transform of this difference equation yields the following transfer function H(z) through the filter:

$\begin{array}{lcl} H(z) & = & \frac{Y(z)}{X(z)} \\ & = & \frac{1}{1 - z^{-1} - z^{-2}} \\ & = & \frac{z^{2}}{z^{2} - z - 1} \end{array}$

This transfer function has two zeros at the origin, and a pair of real-valued poles at z=-0.618 and z=1.618. A pole-zero plot can quickly be produced using the roots and zplane commands in Matlab:

zplane(roots([1,0,0]), roots([1,-1,-1]))


For the transfer function H(z), the coefficient matrices of the numerator and denominator polynomials in z-1 are given by B=[1] and A=[1,-1,-1], respectively. The frequency response of the filter can easily be computed using the freqz command in Matlab. In the following snippet of code, H is the complex-valued DFT and W is a vector of radian frequencies at which H has been computed. The following plot shows the magnitude of the frequency response.

% numerator polynomial
B = [1];
% denominator polynomial
A = [1,-1,-1];
% compute frequency response
[H,W] = freqz(B,A);
% plot the magnitude response (normalize W to pi)
plot(W/pi,abs(H));


The frequency response is computed along the unit circle in the Z-domain (z=ejw). Since the response is symmetric, the above plot only shows W/pi going from 0 to 1. As expected, the magnitude response shows peaks near 0 and 1, due to the relative proximity of the poles to these two points.

In the time-domain (or sample-domain), the response of such a system to a particular input sequence can easily be computed using the filter command in Matlab. In the following fragment of code, the above filter polynomials are used, along with initial states of zero. The input stimulus sequence x is ten samples of a Kronecker delta function, starting at n=0 and ending at n=9. The filter output y is initially 0, but then proceeds through the next ten numbers in the Fibonacci sequence.

% Kronecker delta function starting from n=0 to n=9
x = [1,zeros(1,9)];
% initial conditions of filter
init = [0,0];
% compute the filter response
y = filter(B,A,x,init);


A quicker way to compute the impulse response is through the impz command. The following line of code not only gives the first 10 values of the impulse response, it also produces a stem-line plot of the impulse response.

% compute impulse response
impz(B,A,10)


A simpler way to implement a Fibonacci sequence generator would be to discard the input x and initialize the filter to [1,0] instead of [0,0]. In any case, the Fibonacci filter example is not really practical, as the numbers increase until the internal states exceed the numerical limitations of the implemented machine. An unbounded response indicates that the filter is unstable, and this is confirmed by the location of a pole outside the unit circle in the Z-domain.