

function reverb = make_reverb(beta, b, time, fs, filterBank, fmb)

    noise = 2*rand(1, fs *  time)-1;  % uniform, limited to [-1,1]
    %noise = randn(1, fs *  time); % normal distrb, but unlimited
    noiseb = filterBank_apply( noise, filterBank );
    
%     Dn_e = time * fs / 3; % expected noise density
%     Dnb = calc_densities(noiseb, noiseb, fmb, fs);
%     noiseb = noiseb .* sqrt(Dn_e ./ Dnb); % equalise the noise bands to the expected value

    %figure
    reverb = zeros(1, length(noiseb));
    for n = (1:length(beta))
%        reverb = reverb +  noiseb(n,:);
        reverb = reverb + b(n) * noiseb(n,:) .* exp( -beta(n) * [1:length(noiseb)]/fs );

        %plot(abs(fft(noiseb(n,:))));  hold on;
    end
    
    
    % Equalise reverb to expected equalisation
    Dr_e =  fs /6 .* b .* b ./ beta; % expected reverb density
    reverb_b = filterBank_apply( reverb, filterBank );
    Drb = calc_densities(reverb_b, reverb_b, fmb, fs);
    reverb_b = reverb_b .* sqrt(Dr_e' ./ Drb);
    reverb = sum(reverb_b);
    
    
%     figure
%     plot(Drb)
%     hold on
%     reverb_b = filterBank_apply( reverb, filterBank );
%     Drb = calc_densities(reverb_b, reverb_b, fmb, fs);
%     plot(Drb)
    

    %reverb = noise; % test
    
end




