n = 1024;        % size of the domain on which densities are defined
interval = 16;   % size of the interval on which the unit density is nonzero

u1 = zeros(n,1);              %  u1 is a discretized version of the unit density.
u1(1:interval) = 1/interval;  %  u1 is (1/interval) in positions  1:interval, and otherwise zero.

plot(u1)

% by Matlab convention, the length of the output of conv(a,b)
% is -1 plus the sum of lengths of a and b.  (longer than a or b)

cu2 = conv(u1,u1);
u2 = cu2(1:n);

cu3 = conv(u2,u1);
u3 = cu3(1:n);

cu4 = conv(u2,u2);
u4 = cu4(1:n);

cu5 = conv(u3,u2);
u5 = cu5(1:n);

cu6 = conv(u3,u3);
u6 = cu6(1:n);

sixintervals = 1:(6*interval);
sixunits = sixintervals/interval;

close all
hold on
plot(sixunits,u1(sixintervals),'b')
plot(sixunits,u2(sixintervals),'c')
plot(sixunits,u3(sixintervals),'g')
plot(sixunits,u4(sixintervals),'y')
plot(sixunits,u5(sixintervals),'m')
plot(sixunits,u6(sixintervals),'r')
title('convolutions of the unit density -- using the conv() function')
hold off
pause

%  now do the same thing, but this time with the convolution theorem
u1 = zeros(n,1);
u1(1:interval) = 1/interval;

U = fft(u1);
u2 = real(ifft(U.^2));
u3 = real(ifft(U.^3));
u4 = real(ifft(U.^4));
u5 = real(ifft(U.^5));
u6 = real(ifft(U.^6));

sixintervals = 1:(6*interval);
sixunits = sixintervals/interval;
close all
hold on
plot(sixunits,u1(sixintervals),'b')
plot(sixunits,u2(sixintervals),'c')
plot(sixunits,u3(sixintervals),'g')
plot(sixunits,u4(sixintervals),'y')
plot(sixunits,u5(sixintervals),'m')
plot(sixunits,u6(sixintervals),'r')
title('convolutions of the unit density -- using the Fourier transform')
hold off
pause
