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

@katjav said:

Hey Bassik, thanks so much for your test results.

Apart from the 125 Hz the found RT are comparable indeed. But the magnitude plots of the second case? Should M0002_S01_R02.etm be comparable to Test 2b proc.wav? Even without a reference, the figure should show a similar shape because it is in deciBels, no?

But well, at least it shows that the Pd patch is not complete nonsense. We're on the right track.

I wonder how the magnitude of a full IR is computed in EASERA. For a 20 second chirp the IR it is at least 40 seconds, anticausal and causal part taken together. Or is the analysis limited to the causal part? No, I would think that full IR means the whole thing with the left side included. Pd can not do such long FFT's, 131072 points is the maximum. Probably the magnitude could be computed as an average of several overlapping frames. But the phase information is lost then.

I consider writing an FFT object for Pd that operates on a named buffer, and works with 64 bit floats internally. Maybe it is then possible to do larger FFT's. This is also important for calculating inverse filters (correction filters), since the full IR spectrum is sometimes needed for that. There is also something wrong with [fft~] and [rfft~] phase interpretation: it does not do a centered FFT. So with an IR centered you get these alternating phases instead of a zero-phase representation. Therefore, complex multiplication of spectra will give wrong results. Quite annoying. A new FFT object should have an option to select centered and non-centered FFT so you could work well with centered IR, and with causal part only. Hmm, plans....

Katja

Hello Katja,

Magnitude is expressed in dBFS and for definition our recorded signal should not go above 0 dBFS hence clipping in the analog domain.
So there is something wrong in EASERA when computing our IRs; I haven't figure out yet as the manual does not contain any reference to wav file loading and analysis.

Regarding which part is computed in Easera:
For acoustic measurements (RT, STI etc) only the casual part should be used.
The purpose of the swept sine technique is to get all the reproduction system non-linearities in the anti casual part (leftmost part) of the deconvoluted IR in order to then eliminate them.
So all the analysis should be done on a IR that has the time 0 where the peak of the IR is; also it will be then windowed to eliminate the late part that is not useful to be investigated as it might contain other artifacts of the deconvolution process due to e.g. the noise in the environment during measurements.

I hope I have understood your question.

Glad this keeps you thinking about new objects for PD.

I will review your post on the new chirp formula.

Thank you

Bassik

Hello Katja,

please have a look at the attached file!
Magnitude plots and waterfall plots with our wav file processed (arrival to IR peak moved and windowing) made in Easera; the magnitude is now about right!!!

Probably Audacity applies some processing to the file when exported after processing (just cutting) that lead EASERA to misleading the IR characteristics.

You can see that our IR has its peak normalized to 0 dBFS, did you apply already all our thoughts on IR scaling??

Still need to find out how to import mic calibration in EASERA....

Bassik

http://www.pdpatchrepo.info/hurleur/RT.xls

Hi Bassik,

These comparisons look quite promising! And this is done with the old wiggly chirp. The low frequency ripples over which I'm cracking my brains hardly play a role because the system's response starts at around 100 Hz. I'm curious, what is actually the system under test?

The promising results are of course no reason to stop the quest for the ideal chirp. I have done some fade-in window tests, and indeed they can completely remove the low frequency ripple. But the fade-in window has to be crazy long before it has good result, stretching fully over the first octave. The ripples actually result from the chirp's sudden start. A regular smooth fade-in like the half Hann window I am using now does not effectively counteract this effect. So now I'm thinking of a critically damped oscillator, applied time-reversed as a window. Or maybe a windowed sinc.

About normalization I had some more thoughts. It occurred to me that the direct response will (in most cases) be spreaded over two neighbouring indexes. Ideally there should be just one index representing the indentity diagonal in a Toeplitz convolution matrix. With the two largest peaks in the center of an FFT frame, the group delay can be estimated over a region which has linear phase. Then, the complete IR phase spectrum should be rotated, to center the IR exactly on one index. After this, peak normalization of the IR should be done. Then, the maximum spectrum magnitude must be found. The multiplicative inverse of this maximum provides the normalization factor to prevent clipping. To do all these things, we really need that FFT-XXL object yet to write.

For acoustic measurements (RT, STI etc) only the casual part should be used. The purpose of the swept sine technique is to get all the reproduction system non-linearities in the anti casual part (leftmost part) of the deconvoluted IR in order to then eliminate them.

This is clear. So 'full IR' actually means the causal part for this purpose. I can't remember having heard about STI. Found this DIY article on the topic:

http://www.synaudcon.com/site/articles/a-do-it-yourselfers-guide-to-computing-the-speech-transmission-index/

Seems rather complicated. Should we build this in Pd too?

Katja

More on chirp windowing, good news.

Yesterday I reported about a Hann-fade-in window eating so much of the chirp, before producing good result. Actually that need not be a problem. The chirp can be extended with an extra octave, and this octave is then sacrificed to the fade-in. With 11 octaves in the chirp, and the first octave fully fade-in-windowed, a beautiful spectrum results, with a fairly steep transition band climbing up to precisely the start of the second octave (around 21 Hz at sr 44K1) and with only the faintest remainders of a ripple.

1/11 of the chirp is now spent on a fade-in. For example, almost 2 seconds of a 20 seconds chirp. So what? The spectral result is excellent, and the chirp length may be extended with a few seconds to compensate the lost time. I see no disadvantages so far.

Think I'll implement this long Hann-fade-in by default in my chirp generator.

Katja

So I have this object [expochirp~] finished and done some simple IR measurement on my MacBook internal soundcard (out >> in loopback). Earlier tests had already shown a substantial dip around Nyquist but that seems reasonable as the signal must be filtered to not generate aliases in the discrete domain.

But now the surprise, at least for me. A couple of days ago I uttered the assumption that the direct impulse response will be spreaded over two neighbouring samples. After all, when a pulse is received from acoustic or analog source, it would be very coincidental if it falls exactly on a single discrete index. From my first experiments however, I learned that the pulse is not spreaded over two neighbouring samples, but... over a whole bunch of samples. It is a sin(x)/x figure! (edit: sin(pi*x)/(pi*x) actually). Like spectral smearing in frequency domain. Should have known that....

Ok then, time zero lies somewhere inbetween two indexes. And because of the sin(x)/x smearing, there is no way of exactly taking the causal part of the IR. Everywhere you try to cut it near the highest coefficient, it gives a completely different spectrum. This can not be resolved until the impulse response is time-shifted by a fraction of a sample-interval, to position the direct response exactly on a single index. It's funny that Farina did not mention this issue in his texts on the ESS method.

Katja

p.s. I'll post [expochirp~] and patches soon, whenever I find opportunity to compile on different systems.

@katjav said:

Hi Bassik,

These comparisons look quite promising! And this is done with the old wiggly chirp. The low frequency ripples over which I'm cracking my brains hardly play a role because the system's response starts at around 100 Hz. I'm curious, what is actually the system under test?

Katja

Hello Katja,

the system under test is a big exhibition hall in london (http://www.vam.ac.uk/collections/fashion/galleries/40/index.html) measured using a genelec 8030A monitor and an omni directional mic only.

Keeping it simple to understand teh results quicker.

Regarding your exp chirp computation, I am afraid this goes above my current knowledge so I need to study a bit before giving you an answer or even a comment.

But we should bear in mind that exp Sweep sine technique is done for measuring acoustic spaces (hence the storing of non linearities of the reproduction and recording system in the anti casual part of the IR) and to obtain high quality IR for convolution and auralisation use.

In this case I don't really get the point why you are researching the exact sample value of the peak value.
In my experience, analysis of simple IR don't give substantial different results if the time zero of the IR is not exactly on the peak but probably I am missing something here expecially regarding the signal processing.

base of the exponentianl or logaritmic sweep is that the instantaneous frequency vary exponentially with time (see http://www.acoustics.net/objects/pdf/article_farina02.pdf)
Hope this can be of help.

Also Farina described a method for fast convolution for ambiophonics that can maybe be of inspiration for the deconvolution engine (see http://www.acoustics.net/objects/pdf/article_farina04.pdf)

Regarding STI, it is a measure of the speech intelligibility in a given space and takes into account many factors of which th emost important are the S/N ratio and the frequency masking.

Voice is a very complex mathematical beast as it has a characteristics frequency respone that vary with spl and also a frequency modulation mechanism that we use to articulate words. this is obvioulsy different for different languags but has been standardised in a certain way.

I will send to you the ISO standard that clearly explains th emethod of calculation impplemented in any commercial software in teh next days.

Also I would like to discuss possible development after the chirp issues are resolved.

I believe that it can be useful to implent the following in the patch:

• Sound editor to delete the anti casual part of the IR and move the arrival time of the peak to zero

• Room acoustics parameters calculation (I will send you the relative ISO and I will write more info later)

• Multichannel IR recording and processing including an Ambisonics module; IEM in Gratz (http://iem.at/Members/noisternig/bin_ambi) has done a great work on this. Our implementation would be much simpler and would include A-format to B-Format encoding and additional analysis.

Hope everything is well in life in the meantime

speak to you soon

Bassik

ExpoChirpToolbox.

This will be the name of an IR measurement toolchain, based on the [expochirp~] object. The first two patches are built now:

• [expochirp-generator], generates an exponential chirp of excelllent quality, together with inverse chirp. Deconvolution-test and inspection of the chirp's spectral characteristics within the patch. Chirp & inverse chirp can be stored with metadata in a textfile.

• [IRrecorder], loads a chirpfile as created with [expochirp-generator], applies the chirp to the system under test. Deconvolution, and inspection of spectral characteristics of the impulse response. The raw IR can be stored as 32 bit floating point .wav file.

Next component on my program will be [IReditor], where you can load a raw IR as produced with [IRrecorder] and edit it: trim the useful portion, and eventually window the tail if necessary. The editor should apply the IR as a filter in a soundfile player, so you can check the resulting sound. When these patches work satisfactory, we can think of additional tools: [IRinverter] to produce correction filters, [IRanalyser] for extracting acoustics parameters, [IRambisonic] for multichannel IR...

Attached is ExpoChirpToolbox.zip containing [expochirp~] binaries for OSX and Windows (I'll do Linux later) plus patches [expochirp-generator] and [IRrecorder]. Please post comments or suggestions.

edit: can't attach a file at the moment, download from here;

http://www.katjaas.nl/temp/ExpoChirpToolbox.zip

Katja

Thank you very much for this Katjav.

Much appreciated.

Hope that my post yesterday have been of some help.
I was just reading it today and it contains some incongruence.

Anyway, I will start some testing on the patch asap and let you know.

I will also prepare some more sketches for B-format encoding and more info for the analyzer.

Cheers

Bassik

@bassik said:

In this case I don't really get the point why you are researching the exact sample value of the peak value. In my experience, analysis of simple IR don't give substantial different results if the time zero of the IR is not exactly on the peak

Some more experimenting has learned me that it depends on the system under test. For an IR of the MacBook internal soundcard loopback, the problem is manifest. For an IR of my room, that is not the case. The difference is probably this: the soundcard's spectrum is near to perfect and there is only a direct response. The room test response missed a lot of the low spectrum end (also due to lousy 'test equipment'), and the direct response is buried in reflections because of the small room size. It's funny that the room response is much easier to trim than the loopback response.

Searching the internet I found the name for the phenomenon: fractional delay. And of course implementations for fractional delay, which can also be used to get rid of fractional delay. As far as I can see now, we'll only need it for cases where the direct response has a prominent role. For example, when we want to produce inverse filters for correction of spectral deficiencies in test equipment. We're not at that point yet.

Anyway, I did a patch with a model of fractional delay, which clearly illustrates the problem. (See attachment). And if you do a loopback test with [IRrecorder], you'll see what I mean as well.

edit: bah, again impossible to attach file. (Did I exceed my upload quotum?). It's here:

http://www.katjaas.nl/temp/fractional-delay.pd.zip

Katja

Patches IReditor and IRanalyser are added to the ExpoChirpToolbox.

The editor is a basic sample-editor for sample-precise trimming of the IR, windowing the tail(s), and normalising the maximum spectrum magnitude to 0 dBFS (full scale). The analyser presents spectral info of an edited IR, and you can test the IR as FIR filter with test signal, sound file or soundcard input. Later, the analyser will be extended to do parameter extraction, but since this will again require a lot of research, I want to share the patches now. Together with [expochirp-generator] and [IR-recorder], they already form a complete toolchain for generating your own FIR filters from acoustic environments.

Katja

ps .zip file attached this time (many thanks to the people who repaired this forum feature)

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

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!!!