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

Dear bassik:
I read your post and a question arose in me when I read this:

"Sweep signal are mono signal and they are commonly reproduced through mono or multi channel mono systems (at least for room acoustics measurements)."

Question: Is it possible that some 'comb filter'-like effect appears when reproducing a mono signal through multiple loudspeakers for acoustic measurements? To my humble knowledge, early delays from each source might get to the microphone at different times for each frequency due to source-wall-mic distances, thus creating the 'comb-filter' effect.

Besides, your team work with @katjav and @lead is great. I wish you all success.

sumidero

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.
The set-up was something like this:

PC --- [< c=~~~ (mic)
Loudspeaker |--0.3m--|

The record was set to start just at t0 = 0.3m / (sound speed in m/s) and ended just before the sine train ended to avoid the recording of other noises. The frequencies were stepped up instead of making a sine. Then the inverse calculus was done to measure the loudspeaker dB response.
I did this some years after that to make some sound treatment to a very small room (with many many normal modes!) which I use for audio creation. It sort of worked decently.
I didn't know about Pd at that time, so I created a long train of ((sine(n*w0*t)-silence-)xn freq bins) with Octave, used a DAW for the playing-recording and later recovered the dB response function with Octave again.

Could it be a choice/alternative for what you want to do?

Hugs to everyone...

sumidero.

@sumidero said:

Could it be a choice/alternative for what you want to do?

I have not heard of this technique yet. Do you know if it is mathematically described somewhere? How many frequencies, are those sine shots windowed, how is the inverse calculated, etc?

At the moment we even have problems with a technique which is described in detail, and one of the problems is the lack of precision in 32 bit floating point format, hindering us in making long sweeps.

I have now done a different implementation of log sweep calculation, using Farina's formula instead of the Matlab code. With Farina's formula, the imprecision pain is spread over the array instead of concentrating in the high indexes. Still not perfect but better than before, see patch.

We're still at step one! Unbelievable how time-consuming this is. In the 'Hamburg' article from Bassik's links is a method for amplitude inversion in time domain. I do not understand it rightaway, have to study a bit more on it.

http://www.pdpatchrepo.info/hurleur/logsweep4.pd.zip

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

@bassik said:

the function T-t) is basically the expression of the test signal

That's what I missed last night, thanks.

So the amplitude scaling is simply (w2/w1)^(-n/(L-1)) and this curve is multiplied with the time-inversed sweep. This is what I have implemented in Pd today, and it looks like the result you have in Excel. It also sounds like a sweep but of course the low frequencies are so soft you can't hear them. See patch.

However there is also this scaling constant C which must be applied to compensate power loss. This is the part that I can't get done. Note that frequencies w1 and w2 are expressed as dimensionless frequencies within range 0 till pi in this text. With equation (8) for the sweep I get good result, but equation (10) for the scaling constant gives excessively large results, due to arraylength L in the nominator. Can you figure this out?

If we have the scaling factor right we can start building the deconvolution engine. Then, if tests are promising, we should find a way to generate better sweeps. This can be done in a specialized external [logsweep] for example, using doubles instead of floats and writing results into named arrays. But this is only worth the effort if we don't find more barriers on our way.

Katja

http://www.pdpatchrepo.info/hurleur/logsweep&inv.pd.zip

@katjav said:

@sumidero said:

Could it be a choice/alternative for what you want to do?

I have not heard of this technique yet. Do you know if it is mathematically described somewhere? How many frequencies, are those sine shots windowed, how is the inverse calculated, etc?

As far as I can remember, there were no phase calculations involved in my calculi, just relative amplitudes for each frequency bin, in dB. I had set the loudspeakers outdoors, side by side and heading upwards to the sky, placed the mic a meter above them and got the 'reference curve' trusting the PA system response (a pair of close field reference monitors), then did the same inside the room. Response over difference in dB and hands on craftmanship to make bass traps. I know nowadays that I maybe should have done something better, with better equipment, I know it was sort of rough, but at that time I just wanted a room with a flatter response, asap.

Nowadays my job has gotten a little hectic, sorry for the delay in giving you feedback.
Best regards,

sumidero

@bassik said:

---

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.

Yes, I think it was it. It was shown on an old CGA monitor but it looks very similar.

I think that one of the main drawbacks of that method might be the influence of the transients in each sine train. If you don't cut them out they would introduce artifacts in the spectra. At the time I did my home measurements, I checked these transients throughout the frequency spectrum and made an algorithm to slice them out before doing dB calculi. I wander where that Octave code has gone to...

I wish you all a good week.

sumidero

Made some good progress on the topic at last.

Attached patch has a log chirp generator for max 20 seconds, and computes the inverse. As a test, the chirp is convolved with it's own inverse, using [partconv~]. The result should approximate a single pulse of value 1 at time L = chirplength.

With stop frequency set at sr/2, that pulse is 0.997393 and surrounding values are in the order of 1/1000. With lower stop frequency, results quickly drift away from the ideal case. The imprecision artefacts as we perceived them, seem to do relatively little harm.

I think that the mathematical routines, though not perfect, are far more precise than physical test equipment (speakers, microphones) will ever be.

So it is now a matter of including record / load / save routines, to collect some real-life responses.

Katja

http://www.pdpatchrepo.info/hurleur/logsweep&inv-test.pd

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

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

Included some recording-, save- and load-routines to the log sweep deconvolution, and done a first real world test in my room. Dreadful standing waves oh boy.

Status quo patch attached.

Some weighted amplitude normalisation should be done on the sweep response. Now I have done an easy peak normalisation, but this will always lead to IR's having too low power.

Anyway, it is already possible to see the effect of the log sweep method (by the way this works with a linear sweep as well): non-linear distortion and pollution is shifted left of time zero and thus separated from the impulse response proper. Use vertical zoom to see this.

It would be useful to include an elementary wave editor for trimming. I have now done this with Audacity. Hardoff has done a wave editor in Pd, we could use a subset of it and include vertical zoom.

Katja

http://www.pdpatchrepo.info/hurleur/logsweepdeconv-test.pd.zip

@bassik said:

I found this software made in Java:

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

It's great and it does all you need plus more.

Our Pd patch may not grow to such a professional tool, if only because it would require years of developement. But from an educational perspective a Pd implementation may be more interesting. I remember very well how IR measurement was demonstrated by our dsp teacher in Sonology class. We went through the mathematics, and I was excited about how you can expand and compress time dimension with help of a sweep.

With our Pd patch, the whole process can be visualised. At any point in it, you can test how (im)precise your method and result is. You can then decide whether it is good enough for your purpose, or should be improved. You can read Ben Saylor's [partconv~] code to see how fast convolution is done in practice. The benefits of open source. I see good reason to continue our log sweep deconv project.

Katja

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

@bassik said:

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

The sweep itself has RMS 1/sqrt(2) as it is purely sinusoidal. Deconvolved, it translates as a unit pulse which, as a filter, would have neutral amplitude effect. A sweep response however, fluctuates with peaks and dips over the frequency range. When peak-normalised, it has only dips and it's RMS will be lower than 1/sqrt(2), assuming it still has sinusoidal character. This also translates to an IR with low power. When used as a filter, such an IR can reduce overall signal output substantially. In my first tests, by some 10 dB or more. See what I mean?

Therefore I think there should be an option to somehow normalise RMS of the sweep response / impulse response. I have no clue how this is handled in regular analysis software. Maybe in frequency domain, after the IR is trimmed to it's useful length?

Katja

@bassik said:

there is a strange behaviour in PD, if i turn on dsp the patch doesn't allow me to move the horizontal slides

This quirk has popped up a few times recently on Pd forum

http://puredata.hurleur.com/sujet-4893-number-slider-boxes-freeze

The deconvolution patch is in no way optimized. The [Scope~] showing sweeps is running at highest possible refresh rate. The deconvolution subpatch does crazy lots of fft's continuously. These things need switches to turn them off occasionally.

Katja

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

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

@bassik said:

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

Ah, so normalisation of the IR is done in such a way that the resolution of the integer format is maximally employed. That makes sense.

And then for convolution of the anechoic signal, a second normalisation factor should be applied, in such a way to prevent clipping. This is the difficult part. How to compute this factor from an arbitrary IR? I don't see how it could be computed from the IR array itself, since the convolution output magnitude depends on input frequency. But I can imagine a test method: convolve a sweep with the IR and find the maximum output, then scale down to have maximum output 0.999999 when expressed in floats. Is that how it is done?

All taken together, lots of routines to be written in Pd. 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?

Katja

@katjav said:

How to compute this factor from an arbitrary IR? I don't see how it could be computed from the IR array itself, since the convolution output magnitude depends on input frequency.

Oops, this is of course very simple to find out in frequency domain since the trimmed impulse response will easily fit in one FFT frame. This thread is now getting full with my stupid remarks, sorry.

Katja