% prcrev.m
%
% Matlab script implementation of the PRCRev late field reverberator.
%
% by Gary P. Scavone

clear all;

% Create our impulse signal
x = zeros(1, 80000);
x(1) = 1;

T60 = 2.0;
fs = 44100;

% Do allpass sections first
g = 0.7;
M1 = 341;
M2 = 613;

ap1b = [-g zeros(1, M1-1) 1];
ap1a = [1 zeros(1, M1-1) -g];
ap2b = [-g zeros(1, M2-1) 1];
ap2a = [1 zeros(1, M2-1) -g];

% Apply filters in cascade
x1 = filter( ap1b, ap1a, x );
x1 = filter( ap2b, ap2a, x1 );

% Now do feedforward comb sections
combM = [1557 2137];
gm = 10.^(-3*combM / (T60*fs));

fbc1 = [1 zeros(1, combM(1)-1) gm(1)];
fbc2 = [1 zeros(1, combM(2)-1) gm(2)];

% The easiest way to implement the feebdback combs for this example is to
% filter the input vector in one function call (method 1). However, you
% could also perform the calculations iteratively in a loop using the
% filter() function, in which case you need to specify filter states as
% inputs and outputs to the function (method 2).
method = 2; % or 2
if method == 1 % vector computation
  y = filter(1, fbc1, x1 ) + filter(1, fbc2, x1 );
else % do iteratively in loop
  y = zeros(size(x));
  fs1 = zeros(1, combM(1));
  fs2 = zeros(1, combM(2));
  for n = 1:length(x)
    [tmp1, fs1] = filter(1, fbc1, x1(n), fs1);
    [tmp2, fs2] = filter(1, fbc2, x1(n), fs2);
    y(n) = tmp1 + tmp2;
  end
end

% Plot and play
clf
plot( y' );
sound(y, fs)
