原理参考：Natural sounding artificial reverberation

combFilter.m:

function output = combFilter(delay, gain, input)

fs = 48000;

delaySample = int32(delayTime * fs / 1000);

B = [1 zeros(1, delaySample - 1)];

A=[1 zeros(1, delaySample - 2) -gain];

output = filter(B, A, input);

end

calcCombGain.m:

function gain = calcCombGain(reverbTime, delayTime)

gain = power(10, -3 * delayTime / reverbTime)

end

allPassFilter.m:

function output = allPassFilter(delay, gain, input)

fs = 48000;

delaySample = int32(delayTime * fs / 1000);

B = [-gain zeros(1, delaySample - 2) 1];

A=[1 zeros(1, delaySample - 2) -gain];

output = filter(B, A, input);

end

reverb.m:

function output = reverb(combDelayTime, combGain, allPassDelayTime, allPassGain, input)

y = zeros(length(input), 4);

combOut = zeros(length(input), 1);

for i = 1:1:4

y(:, i) = combFilter(combDelayTime(i), combGain(i), input);

combOut = combOut + y(:, i);

end

allPassOut1 = allPassFilter(allPassDelayTime(1), allPassGain(1), combOut);

output = allPassFilter(allPassDelayTime(2), allPassGain(2), allPassOut1);

output = output * 0.25 + input;

end

main.m:

clc

clear

T60 = 2000;

combDelayTime = [29.23 37.67 41.49 44.31];

combGain = calcCombGain(T60, combDelayTime);

allPassDelayTime = [5 1.7];

allPassGain = [0.7 0.7];

%input = [1, zeros(48000-1, 1)];

[input fs] = wavread('test.wav');

y = reverb(combDelayTime, combGain, allPassDelayTime, allPassGain, input);

wavwrite(y, fs, 'reverb.wav');

figure(1)

plot(y)