close all; clear; clc;

%%
% Make the probability density function
x = linspace(-5,5,1e2); % PDF discretization axis
pdf = (x>-4 & x<-1).*(1+0.3.*cos(pi*x)) + (x>0).*exp(-x);
pdf = pdf./sum(pdf); % Normalise pdf
cdf = cumtrapz(pdf+eps); cdf(end) = 1; % Cumulative distribution function
% Generate random points according to pdf
N = 1e5; % Number of points to generate
y = interp1(cdf,x,rand(1,N),'linear');
% Calculate the histogram
bins = linspace(-5,5,round(sqrt(N))+1); % Histogram bin edges
yh = histcounts(y,bins); % Count density of occurrences in each bin
xh = bins(1:end-1) + diff(bins)/2; % Bin centers

%% Plot
% Split into blocks
temp = {};
for j = 1:5
    temp{j} = mean(reshape(y,[2^(j-1) N/(2^(j-1))]),1);
    temp{j} = temp{j}(1:2^12);
end

figure(1); clf;
cm = inferno(5+1);
s = subplot(1,2,2); hold on;
for j = 5:-1:1
    [yh,xh] = histcounts(temp{j},bins);
    xh = xh(1:end-1) + diff(xh)/2;
    yh = yh./max(yh);
    plot(xh,0.9*yh-(j-1),'color',cm(j,:));
    patch(xh,0.9*yh-(j-1),cm(j,:),'facealpha',0.2,'edgecolor','none');
end
grid on; box on;
xlim([-5 5]); ylim([-4.1 1.1]);
axis square;
s.YTickLabel = [];
xlabel('x');

[yh,xh] = histcounts(temp{j},bins);
xh = xh(1:end-1) + diff(xh)/2;
yh = 0.9*yh./max(yh);

f = @(a) mean((yh-interp1(x,a*pdf,xh)).^2);
af = fminsearch(f,1);
plot(x,af.*pdf,'k-');

% Plot coloured noise
N = 1e2;
s = subplot(1,2,1); cla; hold on;
plot([1 N],[0 0],'k--');
plot([1 N],[0 0]-1,'k--');
plot([1 N],[0 0]-2,'k--');
plot(0.2*randnd(-2,1,N));
plot(0.2*randnd(0,1,N)-1);
plot(0.2*randnd(2,1,N)-2);
axis square; box on;
xlim([0 N]); grid on;
ylim([-3 1]);
s.YTickLabel = '';
xlabel('samples');