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

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

Hello,

it depends on what you would like to transmit between PD and live:

Audio or midi??

For creating virtual midi channels you can use MIDI YOKE (http://www.midiox.com/)

in case of virtual audio channels you can use Jack server for windows (http://jackaudio.org/)

or rearoute from reaper (http://www.reaper.fm/)

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

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

@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

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

@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 Katjav,

As promised I have performed some comaprisons between the IR measured by EASERA software and the one done with the PD patch we are developing.

As said, I have used the old chirp computation patch and have obtained useful IRs only for the first 2 measurements.

comparison in RT results are in the attached file showing:

• RT results
• Waterfall plots (where possible to compute)
• Magnitude plots

Interesting is to see that EASERA is not able to compute the dB SPL of the wav file as it doesn't have any reference of the microphone calibration.
I need to investigate more if I can insert this info for the wav file so for the moment the waterfall plots are not comparable (EASERA measurements are in dB SPL while PD measurements in dB FS).

Anyway this is just a start and you can see that the RT measurements are quite similar apart for the 125 Hz results.

Hope this helps

More to come

Cheers

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

Hello Katjav,

I haven't used the latest chirp~test patch but the previous version of the soft with teh old chirp computation.

So it was supposed to work at least in terms of recording the chirp and deconvolve it; but as i said it needs more testing as it was just a one off mesurements set the one i performed.

I am sorry as well I haven't contributed much in teh last time but I am full of stuff at work.

Anyway I will post a comparison between the PD measurement and a professional soft measurement in teh same positions soon.

Also I am not able to use your latest chirp~test patch as I am on windows.

thank you for keeping the study alive!!!

Hello Katjav,

hope you are well.

I have performed some test on a real case and I have to say that despite the first 2 measurements worked very well, after that the deconvolution process stopped to work properly.

Attached there is one of the wrong results I have obtained from the soft.

you can clearly hear that instead of the IR you get the inverse filter at the same time where the IR should be.

Is it a user problem?
Is it a multi measurement issue with the software?

I will investigate further before any conclusion.

thanks

Hello Katjav,

I am finally back even if still with some problems as I don't have internet at home and work is manic....

Anyway let's put everything in the right order:

1. [chirp~]: I am not on mac but I can test it on mac at home but it would be good to compile it for all platforms...do you have any reference where I can understand how to do the compiling?
   also when you consider linear sweep, are you considering a different inverse filter which only the time reversal of the test signal?

2. In attachment there is an extract of the CATT acoustics manual on calibration of IR and convolution; It is related to audio convolution for auralisation simulations but the same principle can be applied to our purposes to calibrate IR and convoluted responses.

3. Katjav wrote: Oh yeah, there is something else: for proper analysis, the impulse response of the test set (speaker&mic) should be factored out. Taking close-mic sweep response, invert magnitude spectrum, compensate, etc. Is this regularly done in IR measurement practice? If so, how?

Normally in IR measurement practices with Swept sine, the convoluted result will appear in the software with all the non-linearities on its left hand side. In all common software (as far as I know) the removal of the non linearities is done by hand just by changing the start point of the result IR and by windowing the end part of it....after this you start analyzing all the parameters you want to derive from the IR.
I hope I answered your question.

4. Calibration of IR and resulted convoluted IR
   as far as I know there isn't anything in literature that explains exactly the routine for calibrations but the discussed process in the previous posts will do the job with the important task of storing the gain information in order to be able to calibrate the convoluted anechoic file later.

also please have a look at this paper: http://pcfarina.eng.unipr.it/Public/Papers/238-NordicSound2007.pdf

from page 13 to 24 it explains all the issues that can occurs with sweep sine method.

I will try to start the tests soon and to try to find more info in literature.

Hope this helps

Thank you

Bassik

http://www.pdpatchrepo.info/hurleur/CATT_acoustics_IR_calibration.pdf

Hello Katja,

one other comment:

gain formula should be:

GAIN = (2^15)-1/MAX for a 16-bit file.

or GAIN = (2^23)-1/MAX for a 24-bit file

CHeers

Bassik

Hello Katja,

Fortunately I had a conversation with a colleague today about IR normalisation for another reason and probably we got our answer.

In simulation softwares (Acoustics 3D models simulations), the auralisation module when deriving an IR and then the convolution of the anechoic file, applies a gain at the two stages and stores the info in the file.

In the case of the IR it applies a gain in order to have the maximum peak digital value (0.9999999999) and then consequently the same gain to all the rest of the IR (24-bit integer).

When doing the auralisation it applies another normalization gain to prevent clipping of the auralised file and takes into account that auralisation is made at 16-bit integer.

a simple program for the IR normalisation (converting the IR to PCM) can be:

Variable definitions

FP(N) = floating point impulse response

PCM(N) = linear impulse response

M = sample length

Algorithm

MAX = 0

For N = 1 to M

If abs (FP(N))> Max then MAX = abs(FP(N))

Next N

GAIN = (2^16)/MAX

For N = 1 to M

PCM(N) = FP(N) * GAIN

Next N

Store GAIN

As you can see the gain is arbitrary and depends on the IR (or the recorded sweep response) and there isn't a general rule as far as I know.

IN our case we should do these normalisation of the file and store somewhere (metadata in the wav file) the gain applied that would be useful in future applicaitons.

IN the case of multichannel IRs this can be a bit more complex but one thing at time is better.

Hope this helps

Bassik

Hello Katja,

thank you very much for your answer.

Unfortunately in the next couple of weeks I would not be able to be present on the web due to work commitments and moving house with my wife (no internet at home...).

So I will try to check the forum and have a look at the amplitude normalisation but I am not sure I will have enough time to study it properly.

Sorry about that.

I will get in contact as soon as everything will be sorted.

CHeers

Bassik

Hello Katjav,

The new patch is great, it works well and already does more than I was expected to reach in such a small amount of time.

I would like to do some tests and comparison with a commercial software with the same functionality but I have just realised that the laptop with the software installed is not available until monday; nevermind.

I have two questions on the patch:

1. there is a strange behaviour in PD, if i turn on dsp the patch doesn't allow me to move the horizontal slides or better it doesn't visualize the movements while with dsp off there is no problem. It is a bit of a problem as you don't see the computation progress slider moving.
   Do you know why?
   No error messages and the last PD installed.

2. would you please expand a bit more about weighted amplitude normalisation?
   in the case of IRs normalising the main peak of the file will make all the IR information readable, or probably I am missing the point here.

so regarding the to-do list:

• Multichannel: it is quite important it could be multichannel input (i.e. ambisonic mic or stereo pairs) while it is much less important as multichannel output (1 signal would be common to all speakers you are testing, no need for different signal).

• fft analysis of IR: great idea but would be great to derive also additional acoustic parameters, I will provide maths for all the parameters we need to implement (ISO 3382 part 1 and part2, sorry I need to send the papers privately as I would get in trouble if I share on the web).

Also I would like to add in the future, a-format to b-format encoding (http://pcfarina.eng.unipr.it/Public/B-format/A2B-conversion/A2B.htm).
This is why 4 capsules tetrahedral mics can be used to analyse 3D soundscapes including directional analysis (along 3 axes) and understand issues in room such as room modes directions, flutter echoes and first reflections direction among others.

I am glad you have done a clean-up of the scaling constant calculation as in the previous patch it was a bit confusing to me.

I will have a more in depth play around tomorrow or friday...i am also in the middle of moving house...:(

Thank you

Bassik

Also forgot to mention that looking around for another reason I found this software made in Java:

http://www.hometheatershack.com/roomeq/

It does IR measurements as well.

Cheers

Bassik

Hello Katja,

that is great news!!!

I was looking into the scaling factor today but get stopped by work stuff and I am planning to have a look tonight.

I will also test the patch as well.

If you want I can do a quick real world test here at work as I have all the equipment.

In my first attachment there was a sort of implementation of reproduce and recording of the signal (if I remember well), would this be of help??

Bassik

Hello All,

this is becoming even more exciting!!

@Sumidero:

"Question: Is it possible that some 'comb filter'-like effect appears when reproducing a mono signal through multiple loudspeakers for acoustic measurements?"

Assuming you are in an Airport departure concourse.
You have a big multichannel PA system that should provide voice messages to all the people in the concourse.
You will have a guy speaking to a mic and says something like "Flight xxx will depart from gate yyy".
This is a mono signal reproduced by a multi-channel system.

Comb filtering are likely to happen but this is what you would like to measure to understand if the system would work and what are the downturns on having such a big PA.

Sumidero wrote: "I saw, many years ago, a PC-XP based ISA card for acoustic measurement of loudspeakers that instead of using a sine sweep it sent short trains of frequency fixed sines, opening the record just when the wave train was passing through the microphone."

is this system you are referring to?

http://www.mlssa.com/

IF yes, MLS (Maximum Lenght Signal) measurement system have been the industry standard until the sweep sine method comes out (expecially the log swept sine).

According to my knowledge, It is a very good way of measuring IRs but it has some problems in handling the non-linearities of the reproduction and recording chain.

@Katja

I will have a look to your patch and post some comments later, thank you very much.

I had a read through the 'Hamburg' article and he actually give us the exact function for the inverse swept sine which includes the amplitude inversion (or amplitude reduction, depends how you call it).
The way it is written is a bit confusing, I know, but after a bit of thinking you can see it like that:

assuming your funstion has a development in time (t) that goes between 0 and L-1.
the inverse function including time inversal and amplitude scaling is:

x^-1(t)= L-1-t)*(w2/w1)^(-t/L-1)

so if you want to write as written on physics books where your time variable t goes between 0 and T, it is like that:

x^-1(t)= T-t)*(w2/w1)^(-t/T) where 0<t<T

the function T-t) is basically the expression of the test signal (that you find previously in the paper) where instead of n (or t) you got L-1-n (or T-t).

this T-t) is expressed like that:

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

I hope it is much more understandable now; I tried to use the same way of writing it in excel.
If you look at the xls file in the zip file i posted previously the formula is written properly and there is a graph showing the result.

This is very time consuming I know but there is no rush in finishing it as it could be an important tool in the free software community, I reckon.

Have a nice weekend guys

Bassik

Hello Katja,

thanks for the check.

We might be able to live with this issue, all depends on the (de)convolution process result.

If it eliminates the modulations as well we might be fine, if not we would need to find another way of building the sinesweeps.

Would you be able to do the same patch for the inverse signal, please?

I would try to have a go as well whenever I would get a second.

I will never stop to thank you for your interest.

Bassik