Hi All,

first of all, a big thank you for this forum as it is extremely useful expecially for a newbie like me.

I am in the process of creating a patch that could generate an inverse filter for an acoustic test signal (Log Swept sine) and then perform a deconvolution of the recording of the signal in the space.

What I am trying to do i recreating the swept sine from its mathematical formula which give me the exact inverse filter.

the formula is quite complex:

t)= sin(((w1*T)/(ln(w2/w1))*(e^(((T-t)/T)*log(w2/w1))-1))*(w2/w1)^(-t/T)

where t= time variable
T=signal lenght in s
w2= end frequency(normally 20 kHz)
w1= start frequency (normally 20 Hz)

Attached there is an initial patch that doesn't work.

Convolution engine might follow but for the moment I am more interested in generating the inverse filter.

Anyone has had experience with this?

I will definitely share the patch when finished

Thank you

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

Hi Guys,
I was referred to this thread by Serafino Di Rosario, and I will test this PD patch for performing ESS measurements. Everything is very interesting for me, and it seems that Katja did a very good job!

Regarding the problems encountered, I give you these infos:

1. sine-phase-matched sweep. This method is very useful when performing distortion measurements, or computing multiple-order IRS to be used in not-linear convolution processor (for emulating the nonlinearities of a device). For the method to work, it is mandatory that the sine sweep is sine-phased not only at the beginning, but also at the end of each octave. This way, each harmonic-order IR will be phase matched with the linear IR. The provided formulation solves this problem, and it is very good to see it explained here so simply.
   The importance of using a phase-synced exponential sweep was first discovered by Antonin Novak, a Ph.D. student of the universities of Prague and Le Mans.

2. The ripple at low frequencies can be controlled by proper fade-in. The choice of the "optimal" fade law is still a big subject under scientific discussion. Hann windowing is just a very initial, suboptimal approach. I plan to investigate further the choice of the optimal fade law, and publish something on this topic, soon.

3. The concept of cutting away everything before the arrival of the direct sound is wrong, in my opinion. The "silence" before the arrival of the direct sound has a very important physical meaning, it is the "time-of-flight" of the sound, and provides an accurate measurement of the distance between the source and the receiver. Furthermore, it contains "background noise", which is a very important quantity to know, for example when deriving STI from the IR measurement.
   So PLEASE, do not cut away this initial silence! If the IR has to be used as a filter for a convolution-based reverb plugin, the plugin must be intelligent enough for analyzing the IR, and giving the user the possibility to keep this initial silence or cut it away. For example, IR-1 from Waves gives these possibilities. In any case, a measured IR of a room should always contain the time-of-flight... Publishing "pre-cutted" IRs is wrong, and in the long run will cause a lot of troubles...

4. The "Fractional delay", for the same reasons, should NOT be corrected! If the time-of-flight is fractional, good, let's stay with this fact. As pointed out, cutting (time-shifting) improperly the measured IR can alter its spectrum. So, please, keep every measured IR as it comes out from the convolution with the inverse sweep... If the higher-order distortion products are not needed, it make sense to only keep the linear part only, but always starting from the true "zero time". Let's make an example: I generate a 20s-long IR at 48 kHz, that is 960,000 samples.
   The inverse sweep will also be 960,000 samples long.
   I play the sweep, and record the room response for, say, 1,200,000 samples, for being sure of capturing the complete reverberant tail even at the higher frequencies.
   Now I convolve the recorded signal (1,200,000 samples) with the inverse sweep (960,000 samples), and I get a convolved signal which is long 2,159,999 samples.
   If I want to keep a 4s long IR, containing only the linear response, you should throw away the first 959,999 samples, and keep the following 192,000 samples.
   As this signal starts from the true "zero time", the main peak will not be at the very beginning, but delayed of an amount corresponding to the source-receiver distance. If it was 10m, it will be 10/340=0.0294 seconds...

5. For performing efficiently convolution of very long filters (in the example above, the Inverse Sweep was nearly 1 million points) it is advisable to employ a partitioned convolution scheme. That is, the filter is splitted in a number of blocks, so that instead of performing a single, very long FFT , a number of shorter FFTs is performed instead. On my web site you find a couple of papers explaining the partitioned convolution algorithm. This is the same algorithm employed in the well-known BruteFIR open-source program by Anders Torger.

Bye!

Angelo Farina

Hello Angelo,

Welcome to this discussion. It is a joy to hear the opinion and advice of the foremost expert in this field. Our explorations are based on your work. I want to thank you for all the publications which you have so generously published online.

With the info in your post, a major point has now become clear to me: time zero in an IR refers to the start of the measurement, not to the arrival of the (first) response peak. I have been confused about this all the time, and happy to learn that it is actually so simple and logical.

As you may have guessed from the discussion, we learn from trial and error, just as much as from literature and education. I was 'initiated' in the IR measurement topic some years ago by the dsp teacher of the Sonology Department of the Conservatory of The Hague, Peter Pabon. Some of the mathematics were revealed, and I was fascinated. But now that we are trying to build our own Pd implementation in this forum cooperation, we run into the practical details. One of these was the undetermined phase of the chirp end in the regular mathematical formulation, causing excessive frequency ripple in many cases. Unaware of the work of Antonin Novak, I developed the earlier described solution, which now turns out to be a re-invention, like it happens so often.

Regarding the partioned convolution method: fortunately this is implemented in Pd, by Benjamin Saylor. We use this for the deconvolution of the long chirp responses (~ 20 seconds).

Hopefully you'll find time to try ExpoChirpToolbox and post your comments. It is a work in progress and we are working on further extensions.

Katja

Hello Katjav,

I have been chasing the swept sine technique with Farina's Aurora tools on the Cool Edit/Adobe Audition platform. I too found that octave extension of low frequency side only to use as fade in made possible much better result. Also small fade out at the end for non zero end state and all was greatly improved. But applied to getting IR of sound card loopback produces wrong answer. AC coupling of output and input stages produces min phase Butterworth type answer and swept sine with fades always contaminates with its sinc like IR. MLS gives correct answer, but of course brings its own baggage.

Your webpage: http://www.katjaas.nl/expochirp/expochirp.html shows IR result that I've seen on many of my screens while working on the fade issue. Very nicely done indeed.

Farina released plugins for Audacity for convolution and sweep generation. Audacity display doesn't provide zoom tools for display of fine detail. The plugins do provide ready reference.

The solution to the fade issue is to not use them in producing the sweep pair. The answer is in generating at least one of the sweeps in the frequency domain.

Here is link to a great paper covering some of this: "Transfer Function Measurement with Sweeps," SWEN MÜLLER and PAULO MASSARANI I found it here: http://www.melaudia.net/zdoc/comparisonMesure.PDF

Chapter 5 has got the goods.

Turns out that Kirkeby inverse is perfect tool for the job, and many more! To bad it is trapped in Cool Edit/Audition domain. Perhaps Farina will liberate it soon. The perfect place to apply it is before all the trouble begins. An arbitrary broadband signal of 2^n length may be inverted with Kirkeby, and convolution of the pair returns extremely good result. Samples directly next to pulse of 1e-6 easy to obtain. Starting with white noise, results emulating MLS are possible. Starting with exponential sweep a sweep is returned by Kirkeby with circular convolution properties: To get right convolution result, two copies of forward sweep need to be used and cropped to central result. Alternately Kirkeby may be directed to return result 4x longer for direct convolution. I found that a 2x long Kirkeby result produced artifacts. One may also use 2 copies of a 1x Kirkeby to recover IR with single instance of the seed wave.

I communicated this with Farina, and he honored me by placing his confirmations to: http://pcfarina.eng.unipr.it/Public/Wolkoff/

There are six wave files: 2^20 sample length forward sweep and inverse sweep generated without fades; and 1x,2x,4x Kirkeby inverses of the forward sweep. Additionally a time shift Kirkeby is posted. This is a hack I proposed as a partial work around for straight convolution. It is simply cutting first half of 1x Kirkeby and pasting it back at the end. Fading the ends controls bandwidth, and provides kernel for auralization applications.

Cheers!

Hello Barleywater,

Thanks for your comments!!

I work in collaboration with Katja to the patch.

We will definitely review all your suggestions.

One thing to mentions: we are not using a fade-out window for the chirp formulation but we have improved the formula in order to end the chirp at phase zero at Nyquist frequency.

Our tests shows that fade-out windows includes artifacts in the deconvoluted signal.

Anyway, I am planning to create a website/blog for the software and will post detail soon.

Cheers

Hi Barleywater,

Thanks for pointing to that fine MÜLLER and MASSARANI article. The section about creating chirps in frequency domain is particularly interesting. Unfortunately, in Pd we have this show stopper: FFT larger than 2^17 not possible. While we certainly want 20second chirps. And with the obligatory zero-padding, that would be 40second or 2^21 pt FFT.

I've experimented with inverse filtering a lot. For correction of high frequencies defaults in a broadband chirp, a short FFT is sufficient, in my experience. But for a low frequencies inverse filter, a chirp size FFT is really required, because the low frequency coefficients stretch over the whole IR length. Further, it seems to me that one should also correct the inverse chirp, if all deficiencies are to be corrected 'before all the trouble begins', no?

It's not that I'm worried about spectral deficiencies in Expochirp's generator as it is now, but the topic of inverse filtering is fascinating in general. I recently stumbled upon an old article by Alan Oppenheim &c. about inverse filtering of mixed-phase signals. The signal or IR is first split in minimum-phase and maximum-phase components (done in cepstrum). Find the article here if you're interested:

http://www.rle.mit.edu/dspg/documents/AVOHomoorphic75.pdf

Katja

Hello Bassik and Katjav,

Kirkeby inverse function creates inverse chirp from forward chirp, before all the trouble begins....

A long sweep from f1to f2has no more bandwidth than a short sweep from f1to f2. Seemingly small FFT can do a lot of windowing. Here is study of room transfer function using FFT bandwidth limited sweeps:

http://waveformfidelity.com/

The windowing produces super smooth spectrum at expense of computational noise, which of course is well below that of the room, speakers, and microphone. The extra fade in method is very quiet, but has a little ripple.... I had not found this prior to the above investigation.

Barleywater

@Barleywater said:

Kirkeby inverse function creates inverse chirp from forward chirp, before all the trouble begins....

Ah, seems I misunderstood your 'before all the trouble begins'. Thought you were talking about correction of the forward chirp, using Kirkeby inverse filtering. The overshoot in a chirp generated with the traditional formula is detrimental to the crest factor. So, the trouble begins with the forward chirp. As Bassik mentioned, this is in Expochirp solved by using an improved chirp formula. The concept was already formulated in an earlier post on this topic, but if you're interested you can find the mathematical definitions in this paper:

http://www.uni-weimar.de/medien/wiki/PDCON:Conference/Pure_Data_implementation_of_an_ESS-based_impulse_response_acoustic_measurement_tool

Also, the C source code of Pd class [expochirp~] is open, so you can even check how it is implemented if you like.

Thanks Barleywater, for your contribution to the topic, and for pointing to interesting articles. It's good to know that Cool Edit Pro has this FFT size limit as well, and that there are ways to get around it with downsampling /upsampling. Although in Pd this would be less trivial as you need to develop your own interpolation routines.

Katja

you need to click 'post reply'

the 'quick reply' box that comes up by default doesn't let you post patches.

it's annoying, but that's the way the forum software works.

Hello,

You can now check a dedicated blog for ExpoChirp:

http://expochirp.tumblr.com/

We are going to update it with the last versions and all the new feature we are planning to include.

Please send us your comments!!!