Deconvolution with Matlab

[el...@reuna.cl]

Hello,

I'm trying to make a deconvolution using the deconv
function of Matlab, but I always get NaN as answers.

As an example, I recorded the impulse response of
a Sony MUR201 digital reverberation unit, and then
convolved it with a dry voice recording. This way I
obtained a voice with the frequency response and
reverbaration characteristics of the MUR201.
Up to this point, everything works fine. However,
if I try to make a deconvolution in order to get
rid of the impulse response, I get a vector with
NaNs (Not-a-Number).

This are the steps that I follow:

1. I load the sound files previously recorded:

» h=wread16('impulse-response.wav');
» x=wread16('voice-recording.wav');

(I take care that the first values of h and x
are not zero)

2. I make the convolution:

» y=conv(x,h);

3. I verify that the result is what I expected:

» plot(y)
ok.

» soundsc(y,44100)
sounds great!

4. The deconvolution:

» [q,r]=deconv(y,h);
after several hours, I get a 1, a 1.032, a 2, a 10e+216
and the rest are all NaNs.

Please help.

Elias Andalaft.

-----== Posted via Deja News, The Leader in Internet Discussion ==-----
http://www.dejanews.com/rg_mkgrp.xp Create Your Own Free Member Forum

[T. Douglas Mast]

In article <6plm5f$b6e$1...@nnrp1.dejanews.com>,
el...@reuna.cl writes:

> after several hours, I get a 1, a 1.032, a 2, a 10e+216
> and the rest are all NaNs.

Sounds like Matlab is dividing by zero somewhere within its
deconvolution algorithm.

> » [q,r]=deconv(y,h);

Instead, you could try to do this by Fourier deconvolution, i.e. compute
a transfer function and then invert it. The basic idea would be something
like:

z=fft(y)./fft(h);
result = real(ifft(z));

This should be faster, but it will probably blow up too. You can
try a few things to stabilize the deconvolution, for instance:

1. Filter the transfer function z before inverting (plot it first
to see if blowups are limited to high freqencies or whatever, so
you know where to filter it):

z=filter.*z; result=real(ifft(z));

2. Add a small number ("regularization parameter") to the
denominator to estimate the transfer function:

z=fft(y)./(fft(h)+epsilon); result=real(ifft(z));

If this works for very small epsilon (near the numerical precision limit)
then you're home free; otherwise, use the smallest possible epsilon that
gives accurate results. You can also combine the ideas from (1) and (2),
and try other tricks.

Hope this helps,

Doug.

--
Doug Mast, Applied Research Laboratory, Penn State

[wes...@media.mit.edu]

I don't know anything about what Sony's impulse response looks like, but a
REAL room impulse response is non-minimum phase and can not be inverted using
Matlab's deconv routine. Like Doug suggested, you need to use FFT's.

When you do the IFFT, however, you'll need to FFTSHIFT the output. The
inverse to a non-minimum phase filter is anti-causal. By shifting the
output, you'll move the zero-lag tap to the middle of the filter. I think
that should do it!

I wouldn't worry too much about prefiltering or the epsilon parameter that
Doug suggested - you *probably* won't have any divide-by-zero errors - like I
said, I think that the only reason that deconv failed is because the impulse
response is non-minimum phase...

-----

Alex Westner .. wes...@media.mit.edu
Audio Stuff .. MIT Media Lab .. 617.253.0312
http://www.media.mit.edu/~westner .. 617.576.2694

Official NonLinear Site: http://www.voicenet.com/~legion

-----

In article <6pnail$14...@r02n01.cac.psu.edu>,

--
-----
Alex Westner .. wes...@media.mit.edu
Audio Stuff .. MIT Media Lab .. 617.253.0312
http

[baba...@gmail.com]

hi everybody
i want to make blind deconvolution of a sound file by matlab; can you help me?

[Bjorn Gustavsson]

baba...@gmail.com wrote in message <f5873ff4-2e03-49fd...@googlegroups.com>...

> hi everybody
> i want to make blind deconvolution of a sound file by matlab; can you help me?
>

Yes I can help! May I suggest:

help deconvblind
HTH