Yop,

this thesis (http://www.tesisenred.net/bitstream/handle/10803/7264/TAPC.pdf?sequence=1) explains a method for reconstructing a body impulse response (BIR) characterizing the passage from a piezo pickup on a violin to the sound recorded from a mic in an anechoic chamber of the same violin. A slow and long glissando is recorded in parallel from pickup and microphone, and the BIR is calculated. The violin pickup sound can then been convoluted to render a more realistic sound, simulating the miking. Hear samples here : http://www.dtic.upf.edu/~aperez/JASAradiation/

I already started trying to make it with Pd, but I'm not sure it will be possible (besides maximum buffer length discussed in topics like http://puredata.hurleur.com/sujet-4189-swept-sine-deconvolution ). Once the real and complex parts from [rfft~] are stored in buffers, I guess it will be possible to make the weighted sum and other calculations. Should be easier to do that in C or something! Another concern is that I'm not sure to fully understand what does mean (page 96(114)) "An important issue for computing the average is that, since we are working in a bin-by-bin estimation, we must take into account the effects of the window applied to the signals. This fact led us to weight each bin estimation by the energy contribution of the corresponding bin in the input signal spectrum, so that frames for which the energy is low for the bin under estimation (more likely to be affected by window sidelobes) give less influence to the average".
I think that implies to define a relative energy seen from the frame point of view, am I right ? If it's not the case, I'd say that (if a and b are the real and imaginary parts of [rfft~]) magnitude would be the square root of energy (sqrt(a²+b²)), but this energy is not seen as a relative value but absolute.
What do you think of it all?

Of course this topic is not only related to violin, as this procedure is adapted from other works with guitars. See This thesis (http://www.google.be/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CC0QFjAA&url=http%3A%2F%2Fmtg.upf.edu%2Fsystem%2Ffiles%2Fpublications%2FBuccci-Andres-Master-thesis-2011_0.pdf&ei=bNBZUvvCBPON0wWyi4D4DA&usg=AFQjCNHwZ0JDaVCOcIOZZvzINL3rBMw5RA&bvm=bv.53899372,d.d2k). I stumbled upon this material following shawnb links in this topic (http://puredata.hurleur.com/sujet-9637-guitar-effects-designing-hardware), thank you again. Look at the videos on http://www.futureguitarnow.com/forum/index.php?topic=1945.0 . These guys use some softs to get impulse responses playing riffs along their beloved guitar hero and convolve their guitar, even with an arduino ! Inspiring. As is this : http://www.weinreichlabs.com/

Sorry for the long post,
Thanx,

Nau

Nau, I have tried inverse filtering / deconvolution in Pd some time ago, and found that most can be done with existing Pd (-extended) objects. Some remarks:

1 - Pd's max FFT frame size is 2^17, above that you get weird results.

2 - Pd's FFT objects don't assume zero phase centered signal. Say you have all spectrum coefficients at phase zero, you don't get a zero phase IR upon IFFT. So for some purposes you may need to rotate the signal or otherwise alternate the spectrum coefficients' phases.

3 - [bsaylor/partconv~] does partitioned convolution. Use it when you want low-latency convolution. You must have a minimum phase kernel to get low latency! See: http://www.katjaas.nl/minimumphase/minimumphase.html.

The serious issues with deconvolution are in the math itself, rather than in limitations of Pd objects. To mention a few: it requires division, which leads to instable results when spectrum coefficients are small or zero amplitude (which they are at the spectrum low and high end), so you need something like 'Kirkeby inverse filter'. Also, the resulting IR may be too big for it's frame and it's tails will wrap around if it was not sufficiently zero padded. I guess such issues are discussed in the thesis.

By the way does Roger Lewis really do real time convolution with an Arduino (standalone)? Did you find out which Arduino board? Did he share his Arduino code?

edit: upon closer reading I see that Arduino is used as controller, however he uses Raspberry Pi and Beagleboard as standalone convolvers. Yeah that's cool.

Katja

Thank you katjav. You're lighting my future journey...

Nau

Forgot to mention one thing. If a window in time domain is problematic because it affects frequencies of a glissando unequally, you may be able to replace that with convolution in frequency domain. In case of a Hann window this is very simple, as the corresponding convolution kernel has only three points: 0.25, 0.5, 0.25.

I use this approach when converting a zero phase spectrum in [partconvEQ10~] to minimum phase. In that case the filter is defined in frequency domain and I don't want to go to time domain just to do the simple window, so it's done by convolution in frequency domain instead. [partconvEQ10~] is here:

http://puredata.hurleur.com/sujet-9515-band-based-low-latency-fir-filtering

The three point convolution is in subpatch [pd compute-minimum-phase-IR] > [pd convolve-spectrum]. If you need to keep phase information, real and imaginary coefficients could be convolved separately. No need for full complex multiplication - the three point kernel itself is zero phase.

Katja

Thank you again !

Hi there,

many many days that I read articles and study Katjav's website, and I think I improved my perception of these concepts. Now I have a magnitude and a phase (difference) spectrum derived from the article cited up.

But as I can only find descriptions of how to compute a minimum-phase kernel from a zero-phase one (and from its magnitude only, reconstructing phases from it), I really feel stupid but I am not able to push forward and generalize the concept to my context (complex zero-phase spectrum).

I feel at the same time close and far from this implementation, and after having delayed this message I decided to send this help message I'm stuck !

Anyone can give me a hint (strangely I could bet who will answer to my whinings )

Thank you,

Nau

@katjav said:

Nau, I have tried inverse filtering / deconvolution in Pd some time ago, and found that most can be done with existing Pd (-extended) objects. Some remarks:

1 - Pd's max FFT frame size is 2^17, above that you get weird results.

2 - Pd's FFT objects don't assume zero phase centered signal. Say you have all spectrum coefficients at phase zero, you don't get a zero phase IR upon IFFT. So for some purposes you may need to rotate the signal or otherwise alternate the spectrum coefficients' phases.

3 - [bsaylor/partconv~] does partitioned convolution. Use it when you want low-latency convolution. You must have a minimum phase kernel to get low latency! See: http://www.katjaas.nl/minimumphase/minimumphase.html.

The serious issues with deconvolution are in the math itself, rather than in limitations of Pd objects. To mention a few: it requires division, which leads to instable results when spectrum coefficients are small or zero amplitude (which they are at the spectrum low and high end), so you need something like 'Kirkeby inverse filter'. Also, the resulting IR may be too big for it's frame and it's tails will wrap around if it was not sufficiently zero padded. I guess such issues are discussed in the thesis.

By the way does Roger Lewis really do real time convolution with an Arduino (standalone)? Did you find out which Arduino board? Did he share his Arduino code?

edit: upon closer reading I see that Arduino is used as controller, however he uses Raspberry Pi and Beagleboard as standalone convolvers. Yeah that's cool.

Katja

Hi Katja,

I have resumed this project after a summer hiatus and have spent the past 3 days reading up on developments in the Development board world in the past 6 months.
The Arduino Uno which I have been using as a controller interfaced with PD running on my I mac does work in real time using the convolve object.

I am currently deciding whether to run the PD on Beagle Board Cubiboard or Wanda board all of which can host a full linux system. Part of me has been wondering why I am bothering when I discovered that the SM pro Audio V machine can run VST´s in standalone mode on its Linux bases VFX library software. The reason I decided its worth continuing is that whilst the V Machine is midi controlable and a stand alone multi effects system customised with faders and Knobs to control parameters in the VST's a neat little self contained stomp effects box with the Arduino controlling the development board is still compelling and user customisable ( the V Machine is completely customisable too so I am still not 100% convinced that it is sensible but since when should what is sensible get in the way of having so much fun.

Anyway the Open Music Labs Stomp shield kit looks like a nice short cut.

http://wiki.openmusiclabs.com/wiki/StompShield

this also interfaces with the computer in their demo , the MP3 files on the site do the stomp shield more justice.

On my PD code I will post it when I clean it up and make it stable I had a few crashes just before my summer break and will be getting that up and running and cleaned up in due course.

For what I am doing the Bodilizer plug in developed by Ingmar Johanssen back in 2008 is still an absolute marvel, I use it for electric guitar models as well as Acoustic guitar and it is pretty amazing. Zoom have a new commercial pedal out called the A3 which is very nice but not user customisable my objective is to deliver a user customisable box which they can include their own IR's of their own or other peoples instruments. The Arduino Tre looks like its quite interesting coming up in Spring 2014.

http://hackaday.com/2013/10/04/tre-when-arduino-meets-beagle-bone/

Now to read the paper in the OP.

Hi Roger,

thank you for sharing, nice to see you here !

I'm still stuck into my IR calculation patch, but soon I will post a big workbench that could help future discussion/troubleshooting.

Nau

Hi Nau,

Great stuff, look forward to hearing progress.
I have re gained momentum now and hope to make some meaningfull in roads before Christmas.

All the best

Roger

Hi there,

here is my work in progress for this project, a workbench patch.

Basically I record a long single note on my violon, with a piezo pickup and with a gooseneck condenser at the same time, translate to the left the condenser waveform by a few sample to re-sync them. I make a wav out of it, left is piezo and right is condenser. Then I try to implement the statistical approach explained in the thesis mentioned in my first post (see pp97-98). Of course here the input spectrum is very restricted as I only play a note, so let's remember this is a proof-of-concept, anyway. Then I hope I get correct amplitude and phase of a non-causal transfer function that binds my "piezo" input to my condenser "output". As the method presented in the thesis gives phases, I try to naïvely implement a modified version of Katjav's [partconvEQ10~] zero-phase to minphase conversion, using now a complex spectrum rather than a real one. And from then I should be able to reconstruct the minimum-phase IR which is the aim of the whole patch....

I don't even know wheher what I'd like to do is possible or has sense, and I clearly feel like a fool trying to do things without having a full understanding of it. I think I understand the "classical" scheme that makes a minphase kernel out of an amplitude spectrum, but here what I'm trying to do is too tough for me !
Of course errors can occur in my "statistical reconstruction", or in the zero-to-minphase part, or both. I'd like to test the zero-to-minphase part with "well-defined", "non dubious" Z-plane coordinates of a zero-phase spectrum to test it, but I don't see at the moment how to create one.

I am sorry, for sure this workbench patch doesn't work but I'll go on trying and hope someone will have a look at it and give some advice. It is full of graphs that could help see what's going on, and I tried to make it rather clean.

Thank you,

Nau

edit : sorry the uploading on the site seems to lead to corruption, dunno why, so here is a temporary link : http://bitshare.com/files/2ehgag5w/ir_construction_nau.zip.html

http://www.pdpatchrepo.info/hurleur/ir_construction_nau.zip

Nau, did you first make a zero phase (magnitude only) spectrum from your recorded note, before converting to minimum phase? Because, originally it will have a mixed phase spectrum. Converting to magnitudes will distort the frequency response, but if the IR / spectrum is very long this may not harm so much. The distortion affects the lowest frequency bin most, which is probably way below your violin note.

If an impulse response is used for filtering only (in contrast with reverb I mean), you don't need a long kernel in the end. It should be possible to cut it down to 2048 or 4096 samples after all conversions are done.

How do you get an impulse response from a violin, I wonder? This is one aspect I've failed to understand so far, when it comes to measuring the IR from an instrument. If you measure the IR from a room, a signal with flat spectrum (noise or sweep) is applied to that room and the response recorded. But if you play a long note on a violin, how do you go through the full spectrum of the violin? Play one long sweep over all four strings?

By the way I downloaded your .zip file but could not open it, not sure if it is the file or my archive manager.

Katja

Thank you Katjav. I don't fully understand what you're explaining, I'll try to read and learn. Let's try to upload again...

Yeah, mixed phase, dunno why I thought it was minimum in my case...

The thesis proposes to play a long and "forte" glissando on the G string for excitating
the violin, but here I'm just in a "proof of concept" patch that works on a single note, so when it'll work it will be time to extend to a glissando.

Thanx again,

Nau

http://www.pdpatchrepo.info/hurleur/ir_construction_8.pd

Thank you Katjav. I don't fully understand what you're explaining, I'll try to read and learn. Let's try to upload again...

Yeah, mixed phase, dunno why I thought it was minimum in my case...

The thesis proposes to play a long and "forte" glissando on the G string for excitating
the violin, but here I'm just in a "proof of concept" patch that works on a single note, so when it'll work it will be time to extend to a glissando.

Thanx again,

Nau

http://www.pdpatchrepo.info/hurleur/ir_construction_8.pd

Nau I've seen your patch now.

The FFT frame for your glissando should be long, to get good resolution. 2^17 is the longest in Pd. When using multiple frames they must be integrated (sum frames pointwise and calculate mean values). You seem to record a sequence of 2048 pt spectra in a long buffer, that's giving weird result.

About the inversion routine I can't say anything, it's complex and not commented. Anyway if you want to make a short filter from a long IR, make sure that the energy is located around index zero as much as possible. That is, zero phase (magnitudes only). Minimum phase is a 'luxury' to get low latency, and in any case the conversion to minimum phase assumes zero phase as the input.

With a project so complicated you should really go step by step, and test the result of each step before you move on.

Katja

Thank you Katjav,

I'm testing this for days and days, and yes, step by step, as every step is separated from others. But as I said before I wanted to make the more I could before asking for help...

Ok for the long frames, I know that from the beginning, but do you mean that it is not so clever to work with much less points for this "proof of concept" step ? I limited myself to 2048 mainly because I thought it would be wise to stick to your examples and just start from there. The 2048 points frames are calculated and stored later in a long file so that this can be retrieived for statistical calculation. These calculations could be buggy, but I implemented them rather carefully so maybe it is close to correct. Anyway that's more the "any filter spectrum to minphase spectrum" part of the algorithm that I'd like to test ! But as I said before I can't figure out yet how I could feed it with a "well-defined", "non dubious" input spectrum.

I am in the fog when it comes to interpretating the results I get in my tables. I mean that I get amplitude and phase reconstruction, and then I transform it to z-plane coords using [poltocar~], and then try to follow the "steps" described in your webpages and so on. But when you tell me things like

"Anyway if you want to make a short filter from a long IR, make sure that the energy is located around index zero as much as possible. That is, zero phase (magnitudes only). Minimum phase is a 'luxury' to get low latency, and in any case the conversion to minimum phase assumes zero phase as the input."

this is rather frustrating as I can't recognize a mixed-phase from a zero-phase in Pd and of course even less able to find the way to convert this mixed-phase to a zero-phase one ! Sorry for noobness. There's a long way between Kreidler's tutorials and such dsp engineering. I read somewhere (I'm not on the right computer now, in the original thesis or in another one pointing to this one), that they just make the mag/phase algorithm described @pp 97-98, and then make a minphase spectrum out of it, using techniques found in "discrete-time signal processing", by Oppenheim and al., nothing more. So I know this is possible to do it, but I couldn't find the right lines in the right chapter in the pdf I got (I'll buy it, some day). Katjav, your website is full of very interesting explanations, but sometimes it's hard to port it to Pd anyway, you now.

Wouldn't it possible to let the "inversion" part of the patch apart for now and just focus on a "generic" patch that would take a mixed-phase spectrum as inputs and would output a min-phase kernel as output ? Wouldn't this be interesting for many people, beyond my personal "timbre stamping" project ?

I don't ask you for direct debug of my patch, even if the spectral transformations are familiar to you as they are cut/paste/modify snippets of your own code, but if you could point to pages to read (for example to understand how to make the difference between mixed or zero phase under pd), that could be useful to me.

Thank you,
Have a nice day,

Nau

Hi Nau,

Sorry if my comments are not very clear, this is so complicated and abstract stuff and it's hard to put it all to words. Also I did not understand how you use the long arrays and wrongly concluded that you get weird results. Maybe the inversion works as intended. What I mean with 'test every step before you move on', is to leave out the minimum phase conversion for the moment and first check if the inverted spectrum does a good job as a filter.

At the end of [pd compute-raw-inverse] you have a magnitude spectrum and a phase spectrum. If you convert these arrays to real and imaginary part and do ifft, you get the mixed phase IR. If you use the IR as filter kernel with [bsaylor/partconv~] you can check by ear if it does more or less what you expect it to do. Don't bother about latency for the moment. The inversion routine is much more important and interesting than the minimum phase thing.

If, from the audible result, you conclude that the inversion filter is properly implemented, you can go to zero phase and from there to minimum phase. For best results, use a high resolution spectrum (longer fft size than 2048 points). If you set all phases to zero before ifft (meaning you use the magnitudes only), you get a zero phase IR in time domain. Very simple! The zero phase IR is symmetric and has it's energy concentrated around index zero, as you will be able to verify in the plot. Because index zero is not in the middle of the plot you see the IR as two halves, one at the start and one at the end. This can not be used as filter kernel unless you swap the halves, but it would have too long latency to use it anyway.

From the zero phase spectrum (the magnitude spectrum), the minimum phase version can be calculated. Neglect the phase spectrum from the inversion routine. This makes things easier to implement as well. If all goes well and the minimum phase spectrum is ifft'd to time domain, you should see an asymmetric IR with the energy concentrated near index zero. This kernel has a low latency, but only when used with partitioned convolution with a small partition size. If you set partition size in {bsaylor/partconv~] equal to block size, the filter will add no latency at all. For block size 64 this is fairly CPU intensive. The filter kernel should be cut to a minimum acceptable length to make it most efficient.

Katja

@katjav said:

How do you get an impulse response from a violin, I wonder?

I think it's kind of outdated at this point, but Angelo Farina wrote some papers about it.

Thanks acreil.

Here's one such paper:

http://pcfarina.eng.unipr.it/Public/Papers/123-JNMR98.PDF

Outdated indeed regarding the test signal type, but an interesting article nonetheless. A test signal (MLS, but nowadays Farina would definitely use an exponential sweep) was applied to the violin bridge using a force transducer. The article discusses a variety of inversion techniques, not very detailed. It is stated that a mixed phase signal can not be properly inverted, because it contains both the timbre components (magnitudes) and the reverb components (phases).

Timbre components can be inverted but reverb components not. Exponential decay in the IR would become exponential rise in the inverted IR, meaning an instable filter. One more reason to use a zero phase or minimum phase inverse filter instead.

Katja

Yop,

@katjav said:

What I mean with 'test every step before you move on', is to leave out the minimum phase conversion for the moment and first check if the inverted spectrum does a good job as a filter.

At the end of [pd compute-raw-inverse] you have a magnitude spectrum and a phase spectrum. If you convert these arrays to real and imaginary part and do ifft, you get the mixed phase IR. If you use the IR as filter kernel with [bsaylor/partconv~] you can check by ear if it does more or less what you expect it to do. Don't bother about latency for the moment. The inversion routine is much more important and interesting than the minimum phase thing.

Wise, indeed. I wasn't sure at all the mixed phase IR could be used to make a first test, and now that you confirm that I'll try it very soon. I'll have to record glissandi and clean them and resync them to go to the next step, with longer frames of course.

@katjav said:

If you set all phases to zero before ifft (meaning you use the magnitudes only), you get a zero phase IR in time domain. Very simple! The zero phase IR is symmetric and has it's energy concentrated around index zero, as you will be able to verify in the plot. Because index zero is not in the middle of the plot you see the IR as two halves, one at the start and one at the end. This can not be used as filter kernel unless you swap the halves, but it would have too long latency to use it anyway.

It is easy, as your work has showed. Even easy for me ! Are you meaning that setting all phases to zero is the only way to get a zero-phase IR ? I mean, in the thesis they compute phases and say that they perform convolution without say much more (except for the method of convolution itself, p.105). So maybe did they do mixed-phase IR convolution at the end.
Another work, http://www.google.be/url?sa=t&rct=j&q=&esrc=s&source=web&cd=2&ved=0CDkQFjAB&url=http%3A%2F%2Fmtg.upf.edu%2Fsystem%2Ffiles%2Fpublications%2FBuccci-Andres-Master-thesis-2011.pdf&ei=4VyiUpiEI8bChAewmIHQAQ&usg=AFQjCNGQlU_TGyg17jl5W9PEkplDv-o0Tw&bvm=bv.57752919,d.ZG4, page 13, states that they convert the reconstructed amp/phase spectrum to a minphase one, they cite Oppenheim's book as a source, I guess this is the one you used, as they take only the amplitude spectrum as input. So I see no evidence in these texts for using the reconstructed phase to create a minphase IR. Maybe is it nonsense, impossible.

Thank you again, Katjav. I know now what I will be doing the next evenings !

Nau

Funny enough, you posted your last message while I were writing mine. The article you mentioned is interesting, I'll dig into it and its references soon.
Another funny thing: I was to mention Andreas Langhoff (second author of your article)' work because I met him in Cremona more than 10 year ago. He attended the violinmaking school as I did, and both of us are engineers. I already used IR given by some research website (possibly the one of Farina) to convolve with my violin, and it makes a nice spectral shaping, even if this is a rough way of working. I did that from the beginning, when I wrote the first post. Sorry for not mentioning it earlier.
And last funny thing, you seem to answer questions I were only about to write down. Time warpiiiiiing

Thank you,

Nau

Forgot to mention one important thing.

Whenever phases in a signal are shifted, the signal will need more space than it's original frame. For convolution it means, the resultant length is the added lengths of the input arrays minus one sample. If you don't give that space, the signal tails will wrap within the frame and you get circular convolution. To avoid this, zeropadding of the input frames is done: each (windowed) input frame is copied into a larger frame and then processed.

I've been wondering for a long time if zeropadding is also required when an arbitrary signal is converted to a zero phase signal via frequency domain. I think the answer must be 'yes'. It is like convolution with a filter designed to make all phases zero. Therefore, zeropadding is required.

Katja