Swept sine deconvolution
@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
Swept sine deconvolution
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:
-
[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? -
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.
-
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.
- 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
Swept sine deconvolution
@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
Swept sine deconvolution
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:
-
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
-
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
Swept sine deconvolution
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
Swept sine deconvolution
@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
FFT Abstraction
Great. You will need the [tabletool] object if you do not have it already.
http://williambrent.conflations.com/pages/research.html#tabletool
Quadrature compander
Thanks for testing Slur. Your patch shows that qompander~ can not handle excessive DC steps properly. Interesting matter. I have checked that compander~ the not-quadrature version can not handle this as well. It is not the DC level of say 200 dB itself (amplitude 100000 times unity) which is the problem, but the steps, the discontinuities. Although this is not regular use of a compander, the problem points to subtle but irritating artifacts which I could not reproduce in test conditions but which I heard in real life sound processing. I also checked limiter~ against your test conditions. Limiter~ can handle this. Now I am beginning to believe that limiter~'s convolution routine is crucial. After all, such extreme DC steps like in your test patch would be smoothened with convolution.
So, possibly we need limiter~'s convolution combined with qompander~'s curve and amplitude method but compensated for the phase distortion, to get fast and good quality noise reduction / compression / limiting. Can we do a FIR filter with Pd objects?
[clip~ 1e-05 100] can be replaced by [max~ 1e-05] that's right. The minimum value is there to prevent instability with very small numbers division which is more than a theoretical problem. The help patch name is qompander~-help otherwise Pd can not find that patch from the right-click menu.
Katja
Extracting partial information
I've tried FFT and [bp~] and I have encountered the same problems. [bp~] actually came out on top! I'm intrigued by your idea of using [sigmund~] in peaks mode, but must confess I do not have experience using it in such a way. Would you mind demonstrating?
What about [fiddle~]? Or is [sigmund~] more appropriate?
edit: Here's what I have developed using [bp~] so far. I haven't decided how to attenuate a partial just yet. I'm thinking of simply building all the possible configurations, unless there is an easier way to do it. Although, I'm not sure that the Q parameter is cutting it when using [bp~]. Let me know what you think of the patch! Any feedback is most welcome.
Also of interest: http://williambrent.conflations.com/mov/pitchHeight.mov
http://www.pdpatchrepo.info/hurleur/ricky_bp_partial_filter.pd
Pattern recognition for arrays?
Just as a quick pointer...
This problem involves the subjects of correlation and autocorrelation (CC and AC).
To get a correlation you do something like convolution, it's expensive/forceful.
In its very simplest form, say you have a long signal of m samples M, and a shorter pattern n samples long called N that you are looking for. Starting at 0,0 you multiply M[0] x N[0], M[1] x N[1]... M[n] - N[n] and look for a minimum in the first derivative over n samples. That tells you that something close to your test signal N appears in M.
What AC gives you is a measure of the periodic similarity within a window. It's like a convolution of a signal with a time reversed version of itself. If a signal is periodic then the AC is high, if it noisy then the AC is low.
Do some searching on the words "cross correlation" and "auto correlation"