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

Here's an implementation of your formula. Were you looking for this ?

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

Thank you.

I cannot check the patch now as I am not at home.

I will have a look on Monday and let you know

THanks very much for the quick answer.

Cheers

HEllo,

I have not disappeared with your help but I am just too busy to develop this.

Sorry about that, I will post something new as i will have it ready..

Cheers,

@bassik said:

Hi All,
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.

I'm looking into the deconvolution process too (any help from anyone would be great!).

Have downloaded the package from here with this code.
I've tried running it on Octave, hacking around to try and get it to run without any luck.

If the deconvolution was to be implemented in Pd, is it as simple as fft of 'response' (initial delay snipped off) minus fft of 'test sweep'? Or have I got the wrong end of the stick?

@lead said:

If the deconvolution was to be implemented in Pd, is it as simple as fft of 'response' (initial delay snipped off) minus fft of 'test sweep'?

It is slightly more complicated. First invert phase spectrum and magnitude spectrum of original sweep to create deconvolution kernel. Then use [partconv~] from bsaylor library for fast convolution.

Katja

@katjav said:

@lead said:

If the deconvolution was to be implemented in Pd, is it as simple as fft of 'response' (initial delay snipped off) minus fft of 'test sweep'?

It is slightly more complicated. First invert phase spectrum and magnitude spectrum of original sweep to create deconvolution kernel. Then use [partconv~] from bsaylor library for fast convolution.

Katja

Okay, not sure where to start with that -- shot in the dark from a vague MaxMSP memory, maybe process the lists from cartesian to polar coordinates (or visa versa) then *-1 +1?
I ask as I'm trying out things, they all sound a bit wrong so far...

@lead said:

maybe process the lists from cartesian to polar coordinates (or visa versa) then *-1 +1?

Not sure what you mean exactly. Admittedly, my suggestion was rather cryptic as well.

Anyway, I was thinking about the situation of a test signal: how long is that signal? Probably too long to fit in one fft frame of reasonable length. In that case, I don't think you can invert it's phase spectrum in frequency domain. But no problem: just do a time reversal in signal domain, and it's phase spectrum will be inverted automatically. To be explicit here: I mean you must time-reverse the samples of the original test sweep (not the test result signals). Time reversal does not change the magnitude spectrum, only the phase spectrum.

So remains inverting the magnitude spectrum of the time-reversed test sweep. A log sweep does not have a flat magnitude spectrum, so you need to compute the multiplicative inverse of the magnitude coefficients (in frequency domain, using polar coordinates indeed here). Beware of possible zero-magnitude points in the original, they can not be inverted of course. After magnitude inversion, go back to cartesian coordinates and from them to time domain. Now you have your deconvolution kernel in time domain.

The deconvolution kernel is the test sweep with phase and magnitude spectrum inverted. It can be considered a (very long) FIR filter. With this kernel you can do fast convolution of the test results, to deconvolve these. For this, you can use Ben Saylor's [partconv~] object. It does 'partitioned convolution' in frequency domain, with zero padding to avoid circular convolution.

Good luck. Please let us know if you have succes. It is an interesting technique to create impulse responses of acoustic conditions. Expensive reverb simulators probably rely on similar techniques. What is actually your purpose for it?

Katja

Thanks for your help Katja, hopefully I can get there.

My purpose is that I'm moving from a desktop based setup to a laptop after ten years. I've built up quite a lot of software and hardware, that I want to emulate on my new setup and from what i've read via the link I posted previously, this seems like a way to do it.
It would also be really cool to be able to collect reverbs (using a transient method instead) from the various places I get to visit.
I feel a bit late to the party with the idea of impulse responses to be honest, but I'm sure an open source pd implementation is something lots of people would find useful.

Attached is my attempt at understanding what you mean, I rarely work on fft things in Pd so it's probably wrong, but hey I tried! I've included three attempts in the patch, in case I'm close with one of them.

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

Stupid me. I overlooked the fact that magnitude inversion has to be done in one single fft frame as well. For your test signal of 352800 samples that means a 2^19 point fft. Though it need not theoretically be a problem, Pd can not handle fft with that blocksize.

Several other options remain.

1. If you could do with a linear sweep test signal, you don't need to invert magnitude spectrum because a carefully constructed linear sweep has flat spectrum. The time-reversal to invert phases would suffice.

2. Invert the log sweep's spectrum in another software that can do very long fft's.

Will [partconv~] be able to handle long impulse responses? According to the object's c code it can do up to 256 partitions. Your test signal is done in 178 partitions of size 2048 so that should be ok.

What is by the way the formula of your test signal? Did you create it in a Pd patch? I'd like to know.

Katja

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.