# Fibonacci and filter 16

Posted by **Loren Shure**,

Over the years, people have posted many versions of code for calculating Fibonacci numbers. I thought I'd collect a few of the ideas here and start a discussion on the merits of the different versions.

### Contents

### Fibonacci Numbers

Here's a place to get some history about the man and his work. Over the centuries, people have that these numbers seem to come up in the natural world. Try searching for Fibonacci and sunflower; I just got 31,600 hits. Fibonacci numbers are a sequence of numbers that can be derived via a recurrence relationship. Start the sequence with two |1|s. Then, each subsequent number in the sequence is found by adding the two most recent numbers. Mathematically stated:

### Straightforward for-loop

Let's calculate a reasonable chunk of the sequence through the straightforward method suggested by the formula above.

clear nf = 102; tic fibf(1) = 1; fibf(2) = 1; for n = 3:nf fibf(n) = fibf(n-1)+fibf(n-2); end ftimes(1) = toc; % Collect tic-toc times.

And now let's look at the first few numbers to be sure we're getting the expected answers.

fibf(1:10)

ans = 1 1 2 3 5 8 13 21 34 55

### for-loop with Preallocation

You've all heard it before. One reason the former calculation took so much time is at least alleged to be because of growing
the array `fibf` each time through the loop. So now let's use the same loop, but preallocate the output array so it doesn't need to grow
during the calculation. We'll compare all the computation times later on in this post.

tic fibp = ones(1,nf); for n = 3:nf fibp(n) = fibp(n-1) + fibp(n-2); end ftimes(end+1) = toc;

### Recursive Function

The mathematical formula above suggests that we could write a Fibonacci sequence algorithm using a recursive function, i.e., one that calls itself. Here's what I came up with. Note that this version grows an array each time. I might have been able to be clever about this. I doubt the code would be as clear, however. Feel free to post a version to prove me wrong!

`type fibrec`

function f = fibrec(n) %FIBREC Recursive function for n Fibonacci numbers. % Minimize the error checking so we don't bog down in it. I have included it %in comments instead. % if ~isscalar(n) | ~isreal(n) | n<0 | fix(n)~=n % error('ArtBlog:fibrec:MBRealPosInt','N must be a real positive integer') % end if n == 1, f = 1; % First element is 1. return; elseif n == 2 f = [1 1]; % First two elements are 1. else % Call fibrec with previous result and compute next one from it. fm1 = fibrec(n-1); f = [fm1 fm1(end)+fm1(end-1)]; end

And now let's use it to calculate the first chunk of Fibonacci numbers.

tic fibr = fibrec(nf); ftimes(end+1) = toc;

### Use filter to Calculate the Sequence

I have one other idea now, to use the function `filter`. Looking at the help,

`help filter`

FILTER One-dimensional digital filter. Y = FILTER(B,A,X) filters the data in vector X with the filter described by vectors A and B to create the filtered data Y. The filter is a "Direct Form II Transposed" implementation of the standard difference equation: a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) - a(2)*y(n-1) - ... - a(na+1)*y(n-na) If a(1) is not equal to 1, FILTER normalizes the filter coefficients by a(1). FILTER always operates along the first non-singleton dimension, namely dimension 1 for column vectors and non-trivial matrices, and dimension 2 for row vectors. [Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final conditions, Zi and Zf, of the delays. Zi is a vector of length MAX(LENGTH(A),LENGTH(B))-1, or an array with the leading dimension of size MAX(LENGTH(A),LENGTH(B))-1 and with remaining dimensions matching those of X. FILTER(B,A,X,[],DIM) or FILTER(B,A,X,Zi,DIM) operates along the dimension DIM. See also FILTER2 and, in the Signal Processing Toolbox, FILTFILT. Overloaded functions or methods (ones with the same name in other directories) help timeseries/filter.m Reference page in Help browser doc filter

I see that it could be helpful since prior output values within the filtered array can be used to calculate subsequent ones.
I want the latest one to add together the previous two outputs. My input, `x` is the array `[1 1 0 0 ... 0]`. And for each output, I need only consider the current input, assuming I pad my original starting value of 1 with zeros.
Therefore, for this case, `b = 1`. The leading coefficient, multiplying the last output, is 1. And then I want to add the two previous outputs, but in the
equation, these appear with `-a(k)`, so the `a` portion of the filter is `[1 -1 -1]`.

tic x = [1 zeros(1,nf-1)]; a = [1 -1 -1]; b = 1; fibfil = filter(b, a, x); ftimes(end+1) = toc;

Let's be sure I got the right answer.

isequal(fibf, fibp, fibr, fibfil)

ans = 1

### isequal Tip

Did you know `isequal` could take more than two inputs? If you store your items to compare in a cell array, `C`, you can then do this: `passfail = isequal(C{:}):`, taking advantage of the comma-separated list notation.)

### Compare Computation Times

It's time to see how quickly the various methods performed. I have multiplied the times by 1000 to make reading the numbers easier.

fprintf('%s\n\n','Compare Times (in milliseconds)') meths = ['loop ',... 'prealloc ',... 'recursive ',... 'filter ']; fprintf('%s\n\n',meths) fprintf('%2.9f ', ftimes*1000), fprintf('\n')

Compare Times (in milliseconds) loop prealloc recursive filter 0.329650836 0.243606380 2.244419333 0.058666674

### filter is Useful for Vectorization

I've shown one example here in which `filter` was useful for vectorizing code and make it perform well. Another one I know about is code from John D'Errico for calculating a moving window standard deviation. There are probably others. Do you know of other good
uses for `filter`, besides straightforward applications for filtering and convolving signals? Let me know.

Published with MATLAB® 7.2

**Category:**- Efficiency,
- Less Used Functionality

## 16 CommentsOldest to Newest

**1**of 16

Hi,

Other examples of running/moving stats. was provided in CSSM by legendary Peter Acklam:

with the title “Running (moving) sum, mean, variance and standard deviation (Was: Vectorized version of moving or running standard deviation)” on Friday 12 Mar 2004 10:16.

However, for my own applications, where I may find occasional missing data (NaNs), FILTER return all NaNs and I use a for loop (Ouch :-)

/Lars

**2**of 16

Here’s my non-array-growing version of fibrec. I only had to change the two lines in the ‘else’ clause.

Steve

function f = fibrec(n)

%FIBREC Recursive function for n Fibonacci numbers.

% Minimize the error checking so we don’t bog down in it. I have included it

%in comments instead.

% if ~isscalar(n) | ~isreal(n) | n

**3**of 16

Well that didn’t work. I’ll try again but only with the two lines I changed in the ‘else’ clause:

f = [fibrec(n-1) 0];

f(n) = f(n-2) + f(n-1);

Steve

**4**of 16

If you work with Fibonacci numbers (or something similar where recursion might come handy) and have enough memory it might be cunning to do something like this (speeding up things by a factor of much):

function f = fibrec(n)

%FIBREC Recursive function for n Fibonacci numbers.

% Minimize the error checking so we don’t bog down in it. I have included it

%in comments instead.

% if ~isscalar(n) | ~isreal(n) | nlength(F)

if n == 1,

F = 1; % First element is 1.

return;

elseif n == 2

F = [1 1]; % First two elements are 1.

else

% Call fibrec with previous result and compute next one from it.

fm1 = fibrec(n-1);

F = [fm1 fm1(end)+fm1(end-1)];

end

end

f = F(1:n);

**5**of 16

Somethings obviously got lost…

Before ‘if n==1’ there should be two lines with:

persistent F

if n>length(F)

To make above work…

**6**of 16

Your filter code tests fine, but the text seems inconsistent with the code. You write:

> My input, x is the array [1 1 0 0 … 0]. And for each

> output, I need only consider the current input, assuming

> I pad my original starting values of [1 1] with zeros.

In the code, x has only a single 1 followed by zeros. The second one is unnecessary because the filter picks it up from y[1], no?

**7**of 16

Michael,

You are correct and I have fixed the text in the blog to reflect that. Thanks for the correction!

–Loren

**8**of 16

If you don’t want all n numbers, the following can give you a single one: Note, this gives f0 = 1, f = 1, f2 = 2, etc…

Also, at around n = 70 (near bitmax) this starts deviating from the results returned by addition based implementations

function fn = fibonnacci(n)

alpha1 = (1 + sqrt(5)) / 2;

alpha2 = (1 – sqrt(5)) / 2;

c1 = (1 – alpha2) / (alpha1 – alpha2);

c2 = (alpha1 – 1) / (alpha1 – alpha2);

fn = round(c1 * alpha1^n + c2 * alpha2^n);

% in fact the alpha2 term is really small, so you don’t

% need it, except for n = 0;

fn = round(c1 * alpha1^n);

**9**of 16

On my PC, I get the following timings for nf=102 :

prealloc :0.015, filter:0.053ms.

For nf=10000, prealloc:3.12,filter:4.61 ms

That is, the prealloc is faster. I have checked this several times. Can this depend on memory, PC internals and so on?

I found another method described in Roman Maeder’s Mathematica book.

Let m=[1 1;1 0];

The matrix, m^(n-1) gives [f(n) f(n-1);f(n-1) f(n-2)]

where f(n) is fibonacci number for n, f(n-1) for (n-1) and so on.

fib=m^(n-1); fib(1,1) gives the fibonacci number for n. As expected,this computation is fast.

It would be interesting to see how the code can be vectorised with the above method to get the entire series upto n. Any volunteers?

**10**of 16

My prealloc results might have been higher because of the way I ran the code. Did you run the preallocation in a script, function, or directly from the command line? That can make a huge difference sometimes, especially command-line vs. file.

**11**of 16

What’s wrong with using the closed form expression? Even if you want to calculate the whole series of numbers it will be much faster. For example,

tic;

a=(1+sqrt(5))/2;

nf=102;

b=1-a;

F=(a.^[1:nf]-b.^[1:nf])/sqrt(5);

toc;

gives me report

Elapsed time is 0.016000 seconds

Eric

**12**of 16

Correction to the previous post.

In fact, the 2nd time I run the code I get

Elapsed time is 0.000000 seconds

I believe the 1st time I ended up with 16ms because of the Hard Disk (or memory) read, while the 2nd time the program was running from cache.

Eric

**13**of 16

hey loren-

i am interested to hear what you or steve think about the idea of matlab supporting gpgpu (using the stream processing hardware of modern video cards to perform vector/matrix operations up to 50 times faster than on the cpu) from matlab. obvious example related to your article: filtering/fft. wired has a great story on recent exciting results in this field:

http://www.wired.com/news/technology/computers/0,72090-0.html?tw=wn_index_16

(i’m not sending to steve because most of the image processing blog is over my head, so i don’t keep up with it — but it seems like that community would also be very interested.)

-erik

**14**of 16

thankyou for the help, it allowed me to cheat!

**15**of 16

hi Loren,

I just wanted to kmow how matlab calls the recursive function when it is repeated twice as in fibonacci series???

**16**of 16

POONKODI-

It just calls it twice. The first call completes (recursion and all) before the second one gets called. This version of recursion is quite wasteful since it doesn’t store and therefore use values that have been previously calculated.

–Loren

## Recent Comments