Matlab : How to represent a real number as binary

Question

Problem : How do I use a continuous map - The Link1: Bernoulli Shift Map to model binary sequence?

Concept : The Dyadic map also called as the Bernoulli Shift map is expressed as x(k+1) = 2x(k) mod 1. In Link2: Symbolic Dynamics, explains that the Bernoulli Map is a continuous map and is used as the Shift Map. This is explained further below.

A numeric trajectory can be symbolized by partitioning into appropriate regions and assigning it with a symbol. A symbolic orbit is obtained by writing down the sequence of symbols corresponding to the successive partition elements visited by the point in its orbit. One can learn much about the dynamics of the system by studying its symbolic orbits. This link also says that the Bernoulli Shift Map is used to represent symbolic dynamics.

Question :

How is the Bernoulli Shift Map used to generate the binary sequence? I tried like this, but this is not what the document in Link2 explains. So, I took the numeric output of the Map and converted to symbols by thresholding in the following way:

x = rand();
 y = mod(2* x,1)  % generate the next value after one iteration

y =

    0.3295 
if y >= 0.5 then s = 1
else s = 0

where 0.5 is the threshold value, called the critical value of the Bernoulli Map.

I need to represent the real number as fractions as explained here on Page 2 of Link2.

Can somebody please show how I can apply the Bernoulli Shift Map to generate symbolized trajectory (also called time series) ?

Please correct me if my understanding is wrong.

How do I convert a real valued numeric time series into symbolized i.e., how do I use the Bernoulli Map to model binary orbit /time series?

Answer 1

You can certainly compute this in real number space, but you risk hitting precision problems (depending on starting point). If you're interested in studying orbits, you may prefer to work in a rational fraction representation. There are more efficient ways to do this, but the following code illustrates one way to compute a series derived from that map. You'll see the period-n definition on page 2 of your Link 2. You should be able to see from this code how you could easily work in real number space as an alternative (in that case, the matlab function rat will recover a rational approximation from your real number).

[EDIT] Now with binary sequence made explicit!

% start at some point on period-n orbit
period = 6;
num = 3;
den = 2^period-1;

% compute for this many steps of the sequence
num_steps = 20;

% for each step
for n = 1:num_steps

    % * 2
    num = num * 2;

    % mod 1
    if num >= den
        num = num - den;
    end

    % simplify rational fraction
    g = gcd(num, den);
    if g > 1
        num = num / g;
        den = den / g;
    end

    % recover 8-bit binary representation
    bits = 8;
    q = 2^bits;
    x = num / den * q;
    b = dec2bin(x, bits);

    % display
    fprintf('%4i / %4i  ==  0.%s\n', num, den, b);

end

Ach... for completeness, here's the real-valued version. Pure mathematicians should look away now.

% start at some point on period-n orbit
period = 6;
num = 3;
den = 2^period-1;

% use floating point approximation
x = num / den;

% compute for this many steps of the sequence
num_steps = 20;

% for each step
for n = 1:num_steps

    % apply map
    x = mod(x*2, 1);

    % display
    [num, den] = rat(x);
    fprintf('%i / %i\n', num, den);

end

And, for extra credit, why is this implementation fast but daft? (HINT: try setting num_steps to 50)...

% matlab vectorised version
period = 6;
num = 3;
den = 2^period-1;
x = zeros(1, num_steps);
x(1) = num / den;
y = filter(1, [1 -2], x);
[a, b] = rat(mod(y, 1));
disp([a' b']);

OK, this is supposed to be an answer, not a question, so let's answer my own questions...

It's fast because it uses Matlab's built-in (and highly optimised) filter function to handle the iteration (that is, in practice, the iteration is done in C rather than in M-script). It's always worth remembering filter in Matlab, I'm constantly surprised by how it can be turned to good use for applications that don't look like filtering problems. filter cannot do conditional processing, however, and does not support modulo arithmetic, so how do we get away with it? Simply because this map has the property that whole periods at the input map to whole periods at the output (because the map operation is multiply by an integer).

It's daft because it very quickly hits the aforementioned precision problems. Set num_steps to 50 and watch it start to get wrong answers. What's happening is the number inside the filter operation is getting to be so large (order 10^14) that the bit we actually care about (the fractional part) is no longer representable in the same double-precision variable.

This last bit is something of a diversion, which has more to do with computation than maths - stick to the first implementation if your interest lies in symbol sequences.

Answer 2

If you only want to deal with rational type of output, you'll first have to convert the starting term of your series into a rational number if it is not. You can do that with:

[N,D] = rat(x0) ;

Once you have a numerator N and a denominator D, it is very easy to calculate the series x(k+1)=mod(2*x(k), 1) , and you don't even need a loop.

for the part 2*x(k), it means all the Numerator(k) will be multiplied by successive power of 2, which can be done by matrix multiplication (or bsxfun for the lover of the function):
so 2*x(k) => in Matlab N.*(2.^(0:n-1)) (N is a scalar, the numerator of x0, n is the number of terms you want to calculate).

The Mod1 operation is also easy to translate to rational number: mod(x,1)=mod(Nx,Dx)/Dx (Nx and Dx being the numerator and denominator of x.

If you do not need to simplify the denominator, you could get all the numerators of the series in one single line:

xn = mod( N.*(2.^(0:n-1).'),D) ;

but for visual comfort, it is sometimes better to simplify, so consider the following function:

function y = dyadic_rat(x0,n)

   [N,D] = rat(x0) ;                   %// get Numerator and Denominator of first element
   xn = mod( N.*(2.^(0:n-1).'),D) ;    %'// calculate all Numerators
   G = gcd( xn , D ) ;                 %// list all "Greatest common divisor"
   y = [xn./G D./G].' ;                %'// output simplified Numerators and Denominators

If I start with the example given in your wiki link (x0=11/24), I get:

>> y = dyadic_rat(11/24,8)
y =
    11    11     5     2     1     2     1     2
    24    12     6     3     3     3     3     3

If I start with the example given by Rattus Ex Machina (x0=3/(2^6-1)), I also get the same result:

>> y = dyadic_rat(3/63,8)
y =
     1     2     4     8    16    11     1     2
    21    21    21    21    21    21    21    21