Deconvolution question

Thanks a lot deadflagblue. It seems like you've filled my roadmap for two monthes (lol). Since my first goal is learning, getting a working system coming second, this will be:

1) Implement Farina's method for exponential sine sweep deconvolution
2) Compare with direct impulse measurement
3) Use Aliki and compare with above
4) Try to equalize a PC speaker system with inverse filtering plus crosstalk cancelation

Some comments on step 1:

My starting points are Agelo Farina's papers 183-AES24, 154-AES110, AES121 Workshop. I got them being an AES member but you may find them on the internet.

Basically, Angelo exposes the maths for exponential sine-sweeps and shows how it separates inharmonic content when the driver/speaker behaves non linearly. The system IR can be recovered by convolving the recorded response with the test signal's inverse filter, being simply itself equalized and time-reversed.

To cut on time, i use the MATLAB program generate_sinesweeps.m (or the result sinesweeps.wav) given by the RealSimple Project Home Page at the excellent CCRMA. Of course, the sine sweep can be done within Max.

http://ccrma.stanford.edu/realsimple/imp_meas/Sine_Sweep_Measurement_Procedure.html

Max plays the file and poke~ the response into a buffer. The buffer is written to a file, processed by Octave/Matlab which writes the result into a second file binded to a Max waveform~.

http://www2.pescadoo.net/wikimages/max/iMac_ess.png

Above is the sine sweep response. Under, the computed IR. The system in use is my iMac speakers and internal mic, obviously not a setup suitable for step 4 but enough for a test. The picture is only of interest when compared to

http://www2.pescadoo.net/wikimages/max/iMac_ess_clipped.png

where i accidentally overloaded the speakers, producing the harmonics peaks. That's where exponential sweeps shine: non-linearities are rejected into negative times. We can just cut them if interested in first-odrer only, or use the Volterra series expansion as descibed in 154-AES110 (i won't do that).

Here are the steps for deconvolving:

- Read sinesweep and response (same length).

- Equalization: As the exponential sine sweep has a 6db/oct emphasis, a 6db/oct roloff filter (simple differentiator) is needed for equalization. In Matlab/Octave, this is filter ([1,-1,0],[0,0], sinesweep). In Max, this would be a onepole~ or a biquad~ with [1, -1, 0, 0, 0].

- Time-reversal: flipud(sinesweep) in Matlab. In max, this would be mxj buf.Op (reverse).

- Zero-padding: as my test signal contains TWO sweeps, simply zero the second half of both signals [Matlab: sinesweep(L/2+1:L) = 0; response(L/2+1:L) = 0;] Should the buffers contain only one sweep, they'd have to be extended to twice their size. This ensures we do a linear convolution, not a periodic one.

- Convolution: ifft(fft(sinesweep).*fft(response)) in Octave/Matlab. The length of the signal (8 seconds here) would imply a Max external like CNMAT's firbank~ or zippernoise's tconvolution~.

That's it.

I don't publish code for now here because it needs to be further validated (it's easy to get outputs, it's less easy to get correct outputs) and the Max code is pretty uninteresting, made of buffer~, poke~, waveform~ ... The tricky part is the Matlab one (although very short). Some refinements in equalization may be necessary to take in account the band-limited characteristic of the system.

But, who knows, a MaxMSP sine sweep deconvo toolbox may be on the way?

Now heading to step2..........