close all;
clear all;

% number of automata
len=12*3;

% number of iterations
iterations=256;

% automata states
state=zeros(iterations,len);
state2=zeros(iterations,len);

% initial state
state(1,len/2)=1;

%for nn=1:len
%    state(1,nn)=rand>0.5;
%end

rtable=[110 30 150 60 24 52 18 26];

% rule
rule=rtable(floor(rand*7+1));

for tt=2:iterations
    % generate rule table
    ruletable=zeros(1,2^3);
    for nn=1:8
        ruletable(nn)=bitget(rule,nn);
    end
    for ii=1+1:len-1
        rule_n=state(tt-1,ii-1)*2^2+state(tt-1,ii)*2^1+state(tt-1,ii+1)*2^0;
        state(tt,ii)=ruletable(rule_n+1);
    end
    if mod(tt,4)==0
        rule=rtable(floor(rand*7+1));
    end
end

image(state.*128); colormap gray;
drawnow;

for tt=1:iterations
    a=state(tt,:);
    state(tt,:)=0;
    b=0;
    for ii=1:len
        if a(ii)==1;
            b=b+1;
            if b==2
                state2(tt,ii)=1;
            elseif b==3
                state2(tt,ii+12)=1;
            end
        end
        if b==4
            break;
        end
    end
end

subscale=[0 2 4 5 7 9 11];

% quantize octaves to scale
for tt=1:iterations
    %if mod(tt,8)==1
    %    nuscale=mod(subscale+floor(rand(1)*12),12);
    %end
    nuscale=subscale;
    for nn=1:len/12-1
        for ii=1:12
            if state2(tt,ii+(nn-1)*12)==1
                for jj=1:length(nuscale)-1
                    if ii~=(nuscale(jj)+1)&&(ii>(nuscale(jj)+1))...
                            &&(ii<(nuscale(jj+1)+1))
                        state2(tt,ii+(nn-1)*12)=0;
                        state2(tt,(nuscale(jj)+1)+(nn-1)*12)=1;
                    end
                end
            end
        end
    end
end

state=state2;

figure;
image(state.*128); colormap gray;
drawnow;

% sound samplerate
samplerate=44100;

% sound length in seconds
seconds=60*2.5;

% ratio of overlap;
%overlap=0.2;
overlap=1.0;

% segment length
seg_len=floor((seconds*samplerate)/iterations);
seg_len_2=floor(seg_len*overlap);

% sound buffer
wave=zeros(1,(iterations+2)*seg_len_2);

% generate window function
%window=blackman(seg_len);
window=ones(1,seg_len);

% generate frequency scale table
scale=2.^(((1:len)-len/2)./12).*440;

t=0;

for nn=1:iterations
    disp(['Synthesizing segment ' num2str(nn) ' out of '...
            num2str(iterations) '...']);
    % synthesize segments
    for ii=1:seg_len
        % sum in the frequencies
        for jj=1:len
            freq=scale(jj);
            % synthesize
            wave((nn-1)*seg_len_2+ii)=wave((nn-1)*seg_len_2+ii)+...
                    window(ii)*(1/len)*state(nn,jj)*...
                        cos(t.*(freq*(2*pi)*1/samplerate));
        end
        t=t+1;
    end
end

wavplay(wave,samplerate);