Study of reverberation equation

Hi all,

Suppose I have following difference equation:

where x(n) is the input audio, y(n) is the output (echoed) audio, D is the echo delay (in samples), and alpha governs the amount of echo fed back.

Now, for a known x(n) as to be my quantized audio signal, I need to study the effect of different values of alpha and D on y(n) i.e. I need to be able to hear different y(n) consequently as a result of variations in D and alpha.

This is how I managed to make a sound on a particular alpha and D:

fs = 8000;
Recorder = audiorecorder(fs, 16, 1);
record(Recorder);
pause(2);
stop(Recorder);
x = double(getaudiodata(Recorder,'int16'));
x = x / max( abs(x) );
b=[1 0 0 0 0 0  0];
a=[1 0 0 0 0 0 -0.4];
y = filter(b,a,x);
sound(y,fs);

But, as mentioned before, I would like to study the impact of different values of alpha and D on y(n) perhaps in something like a for loop. Can anyone assist?

[Guillaume]

The filter you've implemented would be:

-0.4*y(n) = x(n-6) + y(n-6)

Not really the same as your equation.

[Hasan Ghorbani]

You were right. Coefficients were in reverse order. I edited the post. Thanks for pointing it out.

Answer 1

[Guillaume]

If the equation you want to implement is indeed

y(n) = x(n) + alpha*y(n-d)

then the proper filter would have coefficients:

b = 1;
a = [1, zeros(1, d-1), -alpha];

You could implement a generic filter function for this:

filteralphad = @(x, alpha, d) filter(1, [1, zeros(1, d-1), -alpha], x);

and use it any way you want, e.g. in a loop to try different alpha and d:

for alpha = -1:0.1:1
  for delay = 1:6
     y = filteralphad(x, alpha, delay);
     %do something with y
  end
end

[Hasan Ghorbani]

Thanks for the reply but what if D is a decimal value e.g. D=0.2? Then your a vector is no longer valid.

zeros(1, d-1)

[Walter Roberson]

If your delay can be fractional in your sampling frequency then you need to look at Fractional FFT and related concepts.

In the case where the delay can be expressed as a rational fraction of the sampling frequency then the input could be fft'd and the result ifft'd at a higher frequency such that the delay became integral in the new frequency. The interpolated data could then be put through the filter, and then the result down sampled by fft/ifft. There could well be a nice way to implement in the frequency domain so only one fft/ifft was required; I am not familiar with signal processing techniques like that.

But if the delay is fundamentally irrational like π then I do not know what you could do.

Answer 2

[Walter Roberson]

D = 6;
b = [zeros(1,D), 1];

alphavals = -1 : .1 : 1;    %I figure negative corresponds to echo cancellation
num_alpha = length(alphavals);

for K = 1 : num_alpha
  alpha = alphavals(K);
  a = [alpha, zeros(1,D-1), 1];
  y = filter(b,a,x);
  if K ~= 1
     pause(2);    %just to give an audible break, you can change or eliminate this
  end
  fprintf('alpha = %g\n', alpha);
  player = audioplayer(y, fs);
  playblocking(player);
  delete(player);
end

Using audiplayer with blocking is to avoid the overlap in sounds that would otherwise occur in a loop. sound() is asynchronous and does not wait for the playing to finish, so if you loop sound() then you end up with multiple sounds at the same time. audioplayer allows us to block to ensure we only have one sound at a time.