function y = reverb2(x,TR,P,Fe)
%        y = reverb2(x,TR,P,Fe)
% reverberation de x via reflexions, echos et passe-tout
% modele de Moorer
% TR = duree de reverberation
% P = poids du son reverbere par rapport a x
%-------------------------------------------------------

if min(size(x)) > 1
    error('reverb de marche qu''en mono')
end

d = [3 5 7 11 13 17];
r = [51 57 61 67 71 79]/1000;
c = 340;

% premieres reflexions + x
ld = length(d);
g = 1./d;
ret = d/c;
u = x;
for k = 1:ld
    uk = (-1)^k*g(k)*retard(x,ret(k),Fe);
    u = somme(u,P*uk);
end

% echos multiples
g = 10.^(-3*r/(TR+1e-16));
lg = length(g);
z = [];
for k = 1:lg
    %w = echfiltre(u,r(k),g(k),Fe)/lg;
    % variante :
    w = echfiltre(x,r(k),g(k),Fe)/lg;
    z = somme(z,w);
end

% passe-tout et sortie
z = allpass(z,.006,.7,Fe);
y = somme(u,P*z);

% attenuation de la queue (pour enlever le clic)
ly = length(y); df = ceil(ly/5);
y(ly-df+1:ly) = y(ly-df+1:ly).*(1+cos((1:df)*pi/df))/2;
y = y/max(abs(y));


%-------------------------------------------------------
% fonctions internes utilisees par reverb :
%-------------------------------------------------------
function y = retard(x,t,Fe)
% retard de t sec. : y_k = x_{k-m}

x = (x(:))';
m = ceil(t*Fe);
y = [zeros(1,m) x];
%-------------------------------------------------------
function z = somme(x,y);
% addition de 2 vecteurs de longeur differente

lx = length(x);
ly = length(y);
x = (x(:))'; y = (y(:))';
z = zeros(1,max(lx,ly));
z(1:lx) = x;
z(1:ly) = z(1:ly) + y;
%-------------------------------------------------------
function w = echfiltre(u,t,gain,Fe)
% echo multiple sur u, via feedback, retard t
% + filtrage passe-bas

% v_n = u_n + gain*v_{n-m}, w = h*v
% H(z) = 1 / [1 - gain*z^(-m)]

dt = 1/Fe; 
m = ceil(t/dt);
lu = length(u);
v = zeros(1,lu+m);

for  n = m+1:m+lu
    v(n) = u(n-m) + gain*v(n-m);
end

b = fir1(10,.1);
w = filter(b,1,v);
%-------------------------------------------------------
function y = allpass(w,t,gain,Fe)
% passe-tout
% y_n = gain*y_{n-m} + w_{n-m} - gain*w_n
% H(z) = [z^(-m) - gain] / [1 - gain*z^(-m)]
% |H(exp(iwt))| = 1

dt = 1/Fe; 
m = ceil(t/dt);
lw = length(w);
y = zeros(1,lw+m);

for n = 1:m
    y(n) = -gain*w(n);
end
for  n = m+1:lw
    y(n) = gain*y(n-m) + w(n-m) - gain*w(n);
end
for  n = lw+1:lw+m
    y(n) = gain*y(n-m) + w(n-m);
end
%-------------------------------------------------------


