Forums > MaxMSP

Deconvolution question

March 18, 2008 | 3:56 am

Dear all,

I’m slowly trying to work out how to convert a room response (sweep tone recorded in a room) into an impulse response. The consensus seems to be that one convolves the reversed original sweep tone with the recorded sweep tone to produce the IR, but one source has told me this type of convolution requires a DFT. The operation is intended to move all the components of the recorded sweep back to a single point in time, like a Dirac function. Is Max/MSP capable of such an operation? fft~ seems more like an STFT…

Any wisdom on this matter would be much appreciated!


January 9, 2009 | 2:28 pm

Hi demers,

I just came upon your post (I know it’s a long time ago) and I was thinking of doing exactly the same thing. So I was wondering if you got it running in the end?

I have generated the sweeps and their inverse in matlab and I’ve made a quick patch for playback and recording of the sweeps for one source-two receivers (for stereo RIRs) and now I reached the tricky deconvolution point which becomes convolution in our case..
Buffir~ performs linear convolution and is surely not an option with a maximum of 256 coefficients.
Trying to implement it with patcher objects, fft~ breaks the output to real and imaginary part which complicates the thing a bit further but considering that if S is your signal and H is your response then FFT convolution is S(k)*H(k) where the S(k) H(k) are complex DFTs of the s(n) and h(n) in the time domain, BUT padded before the DFT with zeros to equal length (L_s+L_h-1) to overcome the circular convolution problem of the IFFT.
Doing the multiplication with real and imaginary parts you’ve got S_real(k)*H_real(k)+j(S_real(k)*H_img(k)+S_img(k)*H_real(k))-S_img*H_img

This operation can be done with patcher objects but what about the zero padding in real-time??
I would be much easier if there was a non real-time fft object operating on large arrays, I don’t really know anything about real-time operations with blocks.

You think it is doable? I’m really dying to achieve that as it would be a native convolution solution for max to play with and also a deconvolution one for various IR measurements.


January 9, 2009 | 4:00 pm

> I would be much easier if there was a non real-time fft object operating on
> large arrays, I don’t really know anything about real-time operations with
> blocks.
Did you try [jit.fft]?
J-F.

Extreme Time Stretching (update Dec 29th)

http://www.jeanfrancoischarles.com


January 9, 2009 | 4:40 pm

Quote: deadflagblue@gmail.com wrote on Fri, 09 January 2009 07:28
—————————————————-
> This operation can be done with patcher objects but what about the zero padding in real-time??
—————————————————-
Patcher ovelap-add:

max v2;
#N vpatcher 1044 640 1547 999;
#P window setfont "Sans Serif" 9.;
#P comment 235 310 102 196617 — odd frames —;
#P comment 300 296 34 196617 index;
#P comment 231 296 44 196617 window;
#P comment 133 296 34 196617 index;
#P comment 251 181 111 196617 N/2 .. N-1 , 0 .. N/2-1;
#P newex 384 88 27 196617 / 2;
#P message 384 67 49 196617 $1;
#P newex 384 44 48 196617 loadbang;
#P newex 244 239 27 196617 < ~;
#P newex 77 239 27 196617 < ~;
#P newex 244 129 27 196617 +~;
#P newex 244 155 62 196617 %~ $1;
#P newex 77 77 116 196617 count~ 0 $1 1 1;
#N comlet (Signal) Odd index;
#P outlet 308 276 15 0;
#N comlet (Signal) Even index;
#P outlet 141 276 15 0;
#N comlet (Signal) Odd window;
#P outlet 244 276 15 0;
#N comlet (Signal) Even window;
#P outlet 77 276 15 0;
#P comment 85 180 45 196617 0 .. N-1;
#P comment 64 296 44 196617 window;
#P comment 65 309 102 196617 — even frames —;
#P window setfont "Sans Serif" 18.;
#P comment 24 21 325 196626 Overlap-add engine for convolutions;
#P connect 8 0 11 0;
#P connect 11 0 4 0;
#P fasten 15 0 11 1 389 230 99 230;
#P fasten 8 0 6 0 82 217 146 217;
#P fasten 8 0 10 0 82 114 249 114;
#P connect 10 0 9 0;
#P connect 9 0 12 0;
#P connect 12 0 5 0;
#P fasten 15 0 10 1 389 114 266 114;
#P fasten 15 0 12 1 389 230 266 230;
#P fasten 9 0 7 0 249 217 313 217;
#P connect 13 0 14 0;
#P connect 14 0 15 0;
#P pop;

Patcher C*~ (complex multiply):

max v2;
#N vpatcher 962 527 1189 765;
#N comlet I;
#P outlet 113 186 15 0;
#N comlet R;
#P outlet 28 186 15 0;
#P window setfont "Sans Serif" 9.;
#P newex 113 158 44 196617 +~;
#P newex 147 116 27 196617 *~;
#P newex 113 116 27 196617 *~;
#P newex 28 158 42 196617 -~;
#P newex 60 116 27 196617 *~;
#P newex 28 116 27 196617 *~;
#N comlet I2;
#P inlet 164 45 15 0;
#N comlet R2;
#P inlet 130 45 15 0;
#N comlet I1;
#P inlet 60 45 15 0;
#N comlet R1;
#P inlet 28 45 15 0;
#P comment 63 25 85 196617 Produit complexe;
#P connect 1 0 5 0;
#P connect 5 0 7 0;
#P connect 7 0 11 0;
#P connect 3 0 5 1;
#P connect 2 0 6 0;
#P connect 6 0 7 1;
#P connect 4 0 6 1;
#P connect 2 0 8 0;
#P connect 8 0 10 0;
#P connect 10 0 12 0;
#P connect 3 0 8 1;
#P connect 1 0 9 0;
#P connect 9 0 10 1;
#P connect 4 0 9 1;
#P pop;

Patcher convsub~

max v2;
#N vpatcher 833 643 1198 959;
#P origin 0 -2;
#P outlet 272 237 15 0;
#P outlet 51 237 15 0;
#P window setfont "Sans Serif" 9.;
#P window linecount 1;
#P newex 84 186 153 196617 ifft~ $1 $1 $2;
#P newex 178 107 104 196617 fft~ $1 $1 $2;
#P newex 51 107 104 196617 fft~ $1 $1 $2;
#N comlet;
#P inlet 51 44 15 0;
#P inlet 272 44 15 0;
#P newex 84 148 153 196617 C*~;
#P connect 2 0 3 0;
#P connect 5 0 6 0;
#P connect 3 0 0 0;
#P connect 0 0 5 0;
#P connect 3 1 0 1;
#P connect 1 0 4 0;
#P connect 4 0 0 2;
#P connect 4 1 0 3;
#P connect 0 1 5 1;
#P connect 5 2 7 0;
#P pop;

Patcher conv~

max v2;
#N vpatcher 1171 44 1628 387;
#P origin 0 2;
#P window setfont "Sans Serif" 9.;
#P newex 151 123 99 196617 index~ $1;
#P comment 81 71 212 196617 half size : $3;
#P comment 81 59 212 196617 buffer size : $2;
#P newex 67 90 264 196617 overlap-add $2;
#B color 14;
#P outlet 143 274 15 0;
#P newex 143 238 27 196617 +~;
#P newex 50 192 111 196617 convsub~ $2 0;
#B color 14;
#P newex 218 192 111 196617 convsub~ $2 $3;
#B color 14;
#P newex 319 123 99 196617 index~ $1;
#P newex 218 166 27 196617 *~;
#P newex 50 165 27 196617 *~;
#P inlet 23 71 15 0;
#P window setfont "Sans Serif" 18.;
#P comment 67 21 350 196626 fft~ based convolution (size < = 4096);
#P window setfont "Sans Serif" 9.;
#P comment 81 47 212 196617 buffer name: $1;
#P fasten 2 0 3 0 28 153 55 153;
#P connect 3 0 7 0;
#P connect 10 0 3 1;
#P fasten 7 0 8 0 55 223 148 223;
#P connect 8 0 9 0;
#P connect 10 1 13 0;
#P lcolor 7;
#P connect 13 0 7 1;
#P lcolor 7;
#P fasten 6 0 8 1 223 223 165 223;
#P fasten 2 0 4 0 28 153 223 153;
#P connect 4 0 6 0;
#P connect 10 2 4 1;
#P connect 10 3 5 0;
#P lcolor 7;
#P connect 5 0 6 1;
#P lcolor 7;
#P pop;

Suitable for convolutions up to 4096 (2048 impulse response with zero padding). For larger RI sizes, you need to use pfft~


January 9, 2009 | 8:29 pm

Thanx for the quick replies guys!

I haven’t got jitter but I’m thinking of upgrading to max 5+jitter soon (with a student license;-)), jit.fft would be probably ideal for non real-time operations.

Wow Bernier, not one but four convolution patches! I’ll try them as soon as I got home. What are the advantages of using pfft~ over fft~ that allows longer convolutions??

Thanks a ton again


January 9, 2009 | 11:19 pm

Quote: deadflagblue@gmail.com wrote on Fri, 09 January 2009 13:29
—————————————————-
> What are the advantages of using pfft~ over fft~ that allows longer convolutions??
—————————————————-

fft~ has a 4k limit. pfft~ a 64k limit, that is 0.74s at 44.1 for a non circular, zero-padded convolution.


January 10, 2009 | 9:07 am

Dear all

I am very, very interested in this thread, can you share your
findings? I am in the same process now and I have to use commercial
deconvolver which is a pain for sharing the IR after…

thanks

pa


January 12, 2009 | 1:14 am

Bernier thanks for the patch it’s working fine, I tested it with a perfect impulse and I got my impulse response back:-).
Just one question the overlap-add subpatch got me a bit confused, isn’t it there for the zero padding of the input-blocks to length N_h? I thought the overlap-add method is direct (linear) convolution of input blocks then summed together?
I see, so pfft~ can take up to 32768 length IRs? I’ll try to make a version using it.

Tremblap, the deconvolution method depends on the way you are measuring your responses depending on the excitation signal. For the sine-sweep technique it becomes a convolution of the recorded signal with the inverse of the test sweep-signal. You can have these signals pre-generated and turn you laptop to an IR measurement box using maxmsp. The problem is that for large spaces you’ll probably have an IR longer than 0.7sec which is the limit for the above patch using pfft~ subpatchers. So it still has to be performed offline (it’s pretty easy with a numerical software like matlab or octave). Unless you use partitioned convolution which breaks and convolves the impulse response in parts. But I have no idea how it’s done at the moment (the info is pretty scarce) and I don’t know if it can be done with max objects.


January 12, 2009 | 9:36 am


January 14, 2009 | 10:47 pm

Quote: deadflagblue@gmail.com wrote on Sun, 11 January 2009 18:14
—————————————————-
> Just one question the overlap-add subpatch got me a bit confused, isn’t it there for the zero padding of the input-blocks to length N_h? I thought the overlap-add method is direct (linear) convolution of input blocks then summed together?
—————————————————-
The zero-padding/overlap-add algorithm transforms a circular (periodic) convolution into a linear convolution suitable for example to process an infinite real-time signal with an impulse response,as opposed to processing a single buffer ‘one-shot’. The zeropad ensures that the products won’t foldback to lower indices.

The pfft~ version is just a bit different, i’ll copy it at the end of this post.I use the fft~ version for crosstalk cancelation or HRTF processing where IR are short. The pfft~ version can do some nice convoverbs but latency is terrible of course.

Note that 64k is not enough for sinesweep deconvolution, And convolution is only a step in obtaining inverse IR. Most authors use AURORA plugins for CoolEdit (once a free product, now Adobe Audition). Logic and other studio software now have i.r. processing as an integrated part (see Apple’s Impulse Response Utility).

I’ve thinked about a Java extra-long FFT (we are not concerned by efficiency here). I’m playing with Octave, an opensource Matlab-like, and i think it’s a workable choice. Jeremy Bernstein’s shell, and filewatch objects can help integration with Max, Maybe there’s is a Max/Matlab bridge somewhere.

I have tons of links on the subject, let me find them…
http://www.pescadoo.net/wiki/wakka.php?wiki=ConvoLiens

http://www.pescadoo.net/wiki/wakka.php?wiki=ConvoPapers


January 14, 2009 | 10:53 pm

As promised:

Patcher pconvsub~

max v2;
#N vpatcher 777 109 1026 310;
#P origin 62 0;
#P window setfont "Sans Serif" 9.;
#P window linecount 1;
#P newex 47 105 108 196617 C*~;
#P newex 141 57 75 196617 fftin~ 2 square;
#P newex 47 142 108 196617 fftout~ 1 square;
#P newex 47 56 75 196617 fftin~ 1 square;
#P connect 0 0 3 0;
#P connect 3 0 1 0;
#P connect 0 1 3 1;
#P connect 2 0 3 2;
#P connect 2 1 3 3;
#P connect 3 1 1 1;
#P pop;
—————————————
Patcher pconv~

max v2;
#N vpatcher 1170 426 1627 769;
#P origin 0 2;
#P window setfont "Sans Serif" 9.;
#P newex 151 123 99 196617 index~ $1;
#P comment 81 71 212 196617 half size : $3;
#P comment 81 59 212 196617 buffer size : $2;
#P newex 67 90 264 196617 overlap-add $2;
#B color 14;
#P outlet 143 274 15 0;
#P newex 143 238 27 196617 +~;
#P newex 50 192 129 196617 pfft~ pconvsub~ $2 1 0;
#B color 14;
#P newex 218 192 129 196617 pfft~ pconvsub~ $2 1 $3;
#B color 14;
#P newex 319 123 99 196617 index~ $1;
#P newex 218 166 27 196617 *~;
#P newex 50 165 27 196617 *~;
#P inlet 23 71 15 0;
#P window setfont "Sans Serif" 18.;
#P comment 49 21 378 196626 pfft~ based convolution (size < = 65536);
#P window setfont "Sans Serif" 9.;
#P comment 81 47 212 196617 buffer name: $1;
#P fasten 2 0 3 0 28 153 55 153;
#P connect 3 0 7 0;
#P connect 10 0 3 1;
#P fasten 7 0 8 0 55 223 148 223;
#P connect 8 0 9 0;
#P connect 10 3 13 0;
#P lcolor 7;
#P fasten 6 0 8 1 223 223 165 223;
#P fasten 13 0 7 1 156 166 174 166;
#P lcolor 7;
#P fasten 2 0 4 0 28 153 223 153;
#P connect 4 0 6 0;
#P connect 10 2 4 1;
#P connect 10 1 5 0;
#P lcolor 7;
#P fasten 5 0 6 1 324 166 342 166;
#P lcolor 7;
#P pop;
—————————————
Some test patch to play:
(Increase i/o vector size if you hear gaps)

max v2;
#N vpatcher 1170 426 1627 769;
#P origin 0 2;
#P window setfont "Sans Serif" 9.;
#P newex 151 123 99 196617 index~ $1;
#P comment 81 71 212 196617 half size : $3;
#P comment 81 59 212 196617 buffer size : $2;
#P newex 67 90 264 196617 overlap-add $2;
#B color 14;
#P outlet 143 274 15 0;
#P newex 143 238 27 196617 +~;
#P newex 50 192 129 196617 pfft~ pconvsub~ $2 1 0;
#B color 14;
#P newex 218 192 129 196617 pfft~ pconvsub~ $2 1 $3;
#B color 14;
#P newex 319 123 99 196617 index~ $1;
#P newex 218 166 27 196617 *~;
#P newex 50 165 27 196617 *~;
#P inlet 23 71 15 0;
#P window setfont "Sans Serif" 18.;
#P comment 49 21 378 196626 pfft~ based convolution (size < = 65536);
#P window setfont "Sans Serif" 9.;
#P comment 81 47 212 196617 buffer name: $1;
#P fasten 2 0 3 0 28 153 55 153;
#P connect 3 0 7 0;
#P connect 10 0 3 1;
#P fasten 7 0 8 0 55 223 148 223;
#P connect 8 0 9 0;
#P connect 10 3 13 0;
#P lcolor 7;
#P fasten 6 0 8 1 223 223 165 223;
#P fasten 13 0 7 1 156 166 174 166;
#P lcolor 7;
#P fasten 2 0 4 0 28 153 223 153;
#P connect 4 0 6 0;
#P connect 10 2 4 1;
#P connect 10 1 5 0;
#P lcolor 7;
#P fasten 5 0 6 1 324 166 342 166;
#P lcolor 7;
#P pop;


January 15, 2009 | 5:54 pm

I think your pconv~ got pasted twice.


January 16, 2009 | 4:16 am

How stupid!

Well, it was just a "pconv~ Buffer 65536 32768" between adc~ and dac~ as a way to illustrate how to invoke the patcher. #3 should always be half of #2.


January 16, 2009 | 12:50 pm

—————————————————-
Note that 64k is not enough for sinesweep deconvolution, And convolution is only a step in obtaining inverse IR. Most authors use AURORA plugins for CoolEdit (once a free product, now Adobe Audition). Logic and other studio software now have i.r. processing as an integrated part (see Apple’s Impulse Response Utility).
—————————————————-

True, 64K are not enough. There are three convolution externals that I know of that do not have this limitation. One is firbank~ by CNMAT’s suite having an optimised overlap-save conv., the other is tconvolution~ with partitioned convolution and the last is karateconv~ (!!) again with partitioned convolution, I checked quickly the last one with a 2sec IR and I think it sounded quite good.

I’ll propose instead of the non-free audition (shame the excellent aurora plugins are only on audition’s format) either the free ACMus IR acquisition software developed at Sao Paolo University or if you have a linux or OSX machine the excellent Aliki by Fons Andriaensen, another one of the Farina’s team, and one of the Ambisonics experts. Both are doing sinesweep deconvolution and ACMus has also MLS deconvolution.
Tremplap, if you want to try code it yourself here’s some matlab code I found around:
http://www.mathworks.com/matlabcentral/fileexchange/19566
and from ACMUS team:
http://www.mathworks.com/matlabcentral/fileexchange/11392
and for the theory search for farina papers on the subject on google scholar, I think they are the first ones coming up.

——————————–
I’ve thinked about a Java extra-long FFT (we are not concerned by efficiency here). I’m playing with Octave, an opensource Matlab-like, and i think it’s a workable choice. Jeremy Bernstein’s shell, and filewatch objects can help integration with Max, Maybe there’s is a Max/Matlab bridge somewhere.
——————————–

There is a pretty good octave/pd bridge and there is a matlab/max bridge by CNMAT but last time I had a look it was only for MACs and I don’t have one (Max/octave would be a killer combo too..) But maybe that’s the best way to go for a max/pd IR acquisition patch,do the synchronous playback/recording in a patch and buffer the recording into a matlab/octave variable and then run the deconvolution script, all automated into the patch.

——————————–
I use the fft~ version for crosstalk cancelation or HRTF processing where IR are short. The pfft~ version can do some nice convoverbs but latency is terrible of course.
——————————–

Heh glad to hear that, I was doing exactly the same in the university with inverse filtering and crosstalk cancelation but instead of HRTFS I was measuring the actual room response (with the listener inside) with in-ear microphones including some controlled strong reflections. The most interesting outcome was that while if you invert only the HRTFs, strong reflections destroy totally the crosstalk cancelation performance, if you invert the actual room impulse response up to a length of the significant reflections, the performance was getting close to optimal. So measure you binaural RIRs in front of your desktop get the invert filters and you’ve got an excellent virtual surround system! Though in this case the filters you need hsould be much longer (not more than 32768 though.) and a good SNR in your RIRs.

Pheww! that was a long post..


January 16, 2009 | 1:30 pm

> Pheww! that was a long post..

but an amazing one! Aliki is my new friend, will report on this
later this month… (trying to finish a piece right now!)

thanks again

pa


January 19, 2009 | 3:50 am

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……….


January 19, 2009 | 2:21 pm

why build something which already has been built.

why not leave deconvolution to the nineties where it belongs and invent a matter transformer instead.

-110


January 19, 2009 | 5:35 pm

I didn’t know getting impulse responses was so out of fashion!!!

But you can transform matter through deconvolution, no?


February 7, 2009 | 10:32 pm

Hey folks,

I must say I was shocked to see all these responses when I checked back – so long after my original post! I’ll have to trudge through all this now, so I’ll reserve my thoughts for later.

For now, I just wanted to confirm that I did in fact learn Octave and wrote a routine there for producing IRs. Still tweaking it but I think I have more of a handle on it now.

j.y.bernier, am I correct in assuming that your patches are in fact designed as self-enclosed abstractions? I’ve been trying to figure out how to produce those (SDK?), so if you have any tips I will certainly welcome them.

Thank you all for not letting this topic fall into the nether regions…


Viewing 19 posts - 1 through 19 (of 19 total)