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

Lead, I found that Pd can do 2^17 point fft apparently without trouble. If we can construct a log sweep of exactly that size, the inverse can be computed within Pd.

The code and description on ccrma.stanford site in your link imply that deconvolution is done in frequency domain using circular (de)convolution, by complex division with the original sweep spectrum. That is why they use two copies of the sweep in one signal.

We can not do that in Pd because of the fft size limitation. But we can do non-circular convolution with [partconv~]. Therefore we need to create our own sweep and inverse sweep.

The line which I do not understand in the Matlab sweep code is:

t = linspace(0,T-1/fs,fs*T);

I am too lazy to learn Matlab/Octave, so I can't translate this to math or Pd. Does anyone happen to know what it does?

Katja

From Octave help:

-- Built-in Function: linspace (BASE, LIMIT, N)
Return a row vector with N linearly spaced elements between BASE
and LIMIT. If the number of elements is greater than one, then
the BASE and LIMIT are always included in the range. If BASE is
greater than LIMIT, the elements are stored in decreasing order.
If the number of points is not specified, a value of 100 is used.

The `linspace' function always returns a row vector.

For compatibility with MATLAB, return the second argument if fewer
than two values are requested.

Thanks Maelstorm. I should learn Octave.

So linspace creates (fractional) indexes for time t in this case. Well it should be possible to create that log sweep in Pd according to our own specification of length. Plus the complex-inverse via frequency domain. I'll give it a try one of these days when I have time. It's not really my project but I find it interesting stuff.

Katja

Nice one Maelstorm and Katja.

Looking forward to seeing what you come up with Katja, cheers!

Hello all.

It is very interesting reading about people who are trying to make pd work in the same way I am trying to do.

I have left this project apart for a while as I discover the limitation on fft size for pd.

So I would be pleased if We can try to do it together if you agree.

A couple of points:

Log sweep are used in IR measurements, preferred to Lin sweep, as they provide a series of advantages (the best is eliminating reproduction and recording chains non-linearities).

Katjav is totally right when it says that the inverse filter of the log sweep is the time inverse of the test signal with its amplitude reduced by a log law.

In the creation of the signal, linspace function is vital to generate a time vector that control the development of the signal in time.

Formula of the inverse signal (included amplitude reduction) is in the first post if the discussion.

I have all the relative papers if it can be of help (not at moment as I am not at home).

My main issue in pd was to create sweeps longer than 5 sec and their inverse if course.

When back home I can send around what I have done so far and We can discuss it if you like.

Thank you for keeping the conversation alive.

Bassik (Sero)

Hey there Bassik, good to see you back involved, do post whatever you have - as you can see above even my unworking contributions have helped us start to figure out options/routes.

Hello,

thank you for being so welcoming!!

Let's start from the scientific papers:

this is a good overall of all the techniques: http://pcfarina.eng.unipr.it/Public/Papers/238-NordicSound2007.pdf

but the two papers i found much more interesting for our case are:

http://pcfarina.eng.unipr.it/Public/Papers/134-AES00.PDF

(this one contains most of the formulas)
http://dafx09.como.polimi.it/proceedings/papers/paper_26.pdf

The attachment contains all i have done in PD plus the check of the equations in excel.

Hope this helps

Let's keep in touch for developing this as it could a great tool!!

Cheers

http://www.pdpatchrepo.info/hurleur/PD_swept_sine_(PD_forum).rar

Hello Bassik, welcome back to (your own) swept sine deconvolution club, and thanks for links and files.

Today I managed to create a beautiful 2^17 points log sweep in Pd, according to the formula in the docs (patch attached). This is only a first step. Next time the inverse. Together we'll make it, for shure.

Katja

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

Hello Katja,

this is great!!!

thank you very much.

there are still a couple of things to work out by the way:

1. at the end of the sweep there an unusual oscillation (probably my speakers but not sure) and I remember I have read somewhere that in generating the sine sweep there was the need to filter the highest frequency; i will try to find the paper and post it here (hopefully).

2. I tried to change the parameters in your code in order to have a longer sweep (2.1^7 = 6 sec at 44100) but the signal doesn't go till the end but stops earlier; so for sure it is my very newbie atempt to change your code but do you think it would be possible to create longer sweeps in PD?
   It is the main issue in my attempt of coding as normally when measuring acoustic environments you need very long sweep to measure everything you need (about 12 sec sometimes)

THis is great work anyway!!!

I will try to have a more in depth look these days to see if i can help you out with more specific questions.

Thank you

Good Night

I have modified the patch in such a way that start and stop frequency can be set at will, and table length can be 2^17, 2^18 or 2^19. It should be a power of two, to fit exactly in an fft frame.

With longer table length and lower stop frequency, it is easier to hear what actually happens. At the end, the signal stops suddenly of course, and the discontinuity produces a kind of 'blip' effect. Maybe we should just create a very short (64 sample) smooth fade-out there, instead of filtering.

It is now easy to produce log sweeps according to our specifications of frequency and length. But the longer sweeps cannot be handled in frequency domain (for inversion) within Pd. 2^17 really is the largest fft size here.

So, let's see if we can get the whole deconvolution thing done with that size. And if so, we can produce a longer sweep plus inverse in another software. After all you need to compute them only once. The rest of the work is for [partconv~] which handles fast convolution in partitions.

Katja

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

Cripes, this is really coming on! - looking forward to having a dig through when I've got free time.

Hello Katja and lead,

Katja - the new version of the sweep code is great!!!

I have taken some time at work now to calculate the exact exponential of 2 to have 3s, 6s, 12s and 24s sweep and this is included in the attached zip file (I have altered your patch, Katja).

regarding the generated sweep there are still two things to note:

1. at the end you have this 'bleep effect' and I agree with Katja that we can create a very short (64 sample) smooth fade-out

2. by the end of the sweep you get a frequency modulation (oscillation mainly) that is a bit strange; this can be due to the fact that I am listening on headphones the signal and I can have an acoustic "allucination" (typical when listening a mono test signal in a stereo set) or due to something in the code...
   comparing it to the signal at this link it doesn't appear to me

http://www.4horsemen.net/binkster/tracks/track06.zip

So probably we need to have a think about the audio output of the patch.

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

Then regarding the convolution, It would be nice to create a self-working measurement environment in PD (so provide convolution as well) but let's discuss what we can do.

I will investigate a bit more the theory behind the deconvolution of the recorded signal with the inverse of the "pure" text signal and see if we can run the convolution in time domain only instead of occupying the frequency domain.

This is really exciting for me

Thank you very much for your effort in this

Bassik

I am having issues with the forum....

Cannot attach the modified patch, I am afraid.
I would try later or if you would like I can email it.

Ciao

Bassik

Hey Bassik,

I have now better listened to the long sweeps, and found the cause of the modulations. It is, as so often, a 32 bit precision issue. Here are 10 consecutive sample-indexes, as expressed in seconds:

print: 11.3397
print: 11.3397
print: 11.3397
print: 11.3398
print: 11.3398
print: 11.3398
print: 11.3398
print: 11.3398
print: 11.3399
print: 11.3399
print: 11.3399
print: 11.3399

With 64 bit floats and 13 decimal places this would be no problem:

11.33786848073
11.33789115646 etc

Bummer. I don't think we can get around this, with Pd.

Katja

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

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