This analyze plug-in may be used to calculate the decay rate (in dB/sec) of an
exponentially decaying sound.
; Computes decay rate (in dB/sec) by making a least-squares
; fit to linear-to-dB values. Decay rate and RT-60 (time to
; drop 60 dB) are returned as labels.
; First convert sound into rms values. From the Nyquist prompt:
; (force-srate sound-srate (rms s))
;
; Highlight the decay portion–it helps to use Waveform (dB)–
; and then call this routine.
I’ll let the actual programmers dip into that, but I want to know where you found a room with one wall? What kills us is multiple reoccurring echoes, so if you’re after the sterile, stand-alone, scientific experiment, then you’re on the right track.
After you successfully handle analyzing real-world echo decay – without the hand clap, thank-you-very-much – the Masters Thesis or very possibly the Doctoral is the tool to get rid of the echoes.
I always thought you could do a terrific echo cancellation by successive approximation. It might take a while, but eventually you would get there. No matter how you analyze an echo-filled performance, the amplitudes and density of the energy are always – always higher than the performance by itself. Use a wild-ass guess correction and apply it. Did the energy in the performance go up or down? Lather, Rinse, Repeat. Eventually, the software would figure out the size of the room and the distances of the walls enough to do useful correction – by chasing the energy dips.
And yes, as with any cyclical correction like this, you have to know when to stop. You do need the watchdog and loop exit in there somewhere.
So that’s what to work on after you get this analysis thing down.
Congratualations
A nice first Nyquist project.
I think that the three plug-ins could be made into one, can’t they?
Eventually coming up with labels for all octaves.
You’ve taken an really interesting approach with the LPM.
There’s a lot of space for improvement though. There is no need to do the analysis sample by sample.
Firstly, I wonder if the RMS defaults are best suited to the task; how about a seperate control for the measurement interval?
The result will be naturally be in some respect “stepy”. We could maybe use the function “follow” to make the curve smoother (I know, the sound is resampled, but only linearly, which connects the measurement points with straight lines). Or how about using the unresampled RMS and adapt the calculation?
The code would be much more efficient if we worked only with sounds (i.e. the RMS-curves for the various octaves). Maybe it is better explained if I provide the source-code with my comments (2 semikolons):
;; not necessary - I think rj
(defun number-to-string (number)
(format nil "~a" number))
;; no need to re-calculate the length of s
;; "len" already holds length in samples.
;; Thus, (setf nn (truncate (min 1000000 len)))
(setf nn (snd-length s 1000000)) ; integer sample length
;; (float nn) does the same
(setf fn (* 1.0 nn)) ; float sample length
(setf dt (/ 1 *sound-srate*)) ; delta-t
;; Normally, (get-duration 1) returns the dur of s
;; However, the sound could be only 1 milion samples long (see above).
(setf dur (* fn dt)) ; duration in time
;; We could now extract the portion of the sound we need.
;; If we knew the length of the impuls, we could insert this value
;; (setf s (extract-abs imp-end dur s))
;; However, we must take into account that the original
;; impuls has not the same length after RMS.
;; In principle we don't have to initialize the variables.
;; But it may add some readibility.
; Least-squares accumulators. For example, sumy2 is sum of y^2
(setf sumy 0.0)
(setf sumy2 0.0)
(setf sumxy 0.0)
;; We will operate on sounds, hence no loop!
; Iterate through to fill accumulators
(dotimes (ii nn)
;; Simply: (setf val (linear-to-db s))
(setf val (linear-to-db (snd-fetch s))) ; dB values are used
;; we must introduce a new set of functions, which can deal with sounds.
;; an integrator and a fn that gets the last smp.
;; Let's call them "integrator" and "last-smp"
;; (setf sumy (integrator val))
(setf sumy (+ sumy val))
;; (setf sumy2 (integrator (mult sumy sumy)))
(setf sumy2 (+ sumy2 (* val val)))
;; (setf sumxy (integrator (mult val dt (integrator (snd-const 1 0 *sound-srate* dur)))))
(setf sumxy (+ sumxy (* val dt ii)))
)
;; Note: the above functions are calculated such that they result in a sound.
;; But we only need the last value of those and can therefore take them now with
;; (psetq sumy (last-smp sumy) sumy2 (last-smp sumy2) sumxy (last-smp sumxy))
;; However, it is better to take the values directly and thus
;; free the memory from the sounds.
;; I din't examine yet, how the dependencies are, i.e. which values can be taken without
;; taking the whole sound for the calculation.
; Sum of x, and sum of x^2 may be computed directly
(setf sumx (* 0.5 dt (- (* fn fn) fn)))
(setf sumx2 (* (recip 6.0) dt dt (+ (* 2 fn fn fn) (* -3 fn fn) fn)))
; Standard least-squares fit. y = A + Bx. Here B is in dB/sec.
(setf LSD (- (* fn sumx2) (* sumx sumx)) )
(setf LSA (/ (- (* sumx2 sumy) (* sumx sumxy)) LSD) )
(setf LSB (/ (- (* fn sumxy) (* sumx sumy)) LSD) )
; RT60 is time to decay 60 dB.
(setf RT60 (/ -60.0 LSB))
;; the next two lines are much simpler expressed with
;; e.g. (setf slopedur (format nil "Decay: ~a dB/sec" lsb))
(setf slopestr (strcat (string "Decay: ") (number-to-string LSB) (string " dB/sec")))
(setf rt60str (strcat (string "RT-60: ") (number-to-string RT60) (string " sec/60 dB")))
;; Could be expressed somewhat simpler, I guess,
;; especially when labels for each octave should be returned too.
(setf li (cons (list 0 dur slopestr) nil))
(setf li (cons (list 0 dur rt60str) li))
Please bear in mind that I’ve wrote those comments on the fly, so the functions are not tested.
There are the two functions missing:
The function takes as argument also the length of the sound, so nn or fn has to be passed too.
If the values from the integrator are always increasing, we could as well use “peak” to get the last (and highest) sample
You need to be careful using linear-to-db with sounds because dB is a log and the logarithm of numbers <= 0 is not defined, so if the selection being analysed contains absolute silence you will get a “NAN” error (Not A Number).
I asked me the same.
According to sengpiel, both procedures are valuable. If the measurement is sound pressure, RMS will be taken. But sound intensity is also possible. I can’t tell exactly, the page looks like scrambled eggs (for screen readers). http://www.sengpielaudio.com/calculator-RT60.htm
Anyway, that calculation is meant to be for designing rooms.
The following manual tells a little bit about the evaluation of a measurement: http://www.ivie.com/download/IEx5_RT60_Manual.pdf
However, it still doesn’t tell us if the rms or the peak value is taken. Obviously, the program works with 100 samples per second and only 5 values are used to calculate the average slope of the decay.
It is interesting that the RT60 is calculated from 20 and 30 dB intervals, which are extrapolated with the average slope. This makes sense since most rooms have an ambient level of 20 dB and firing pistols to get an 80 to 85 dB impuls might be a bad idea…
I wonder if the evaluation couldn’t be done in the following way:
recording an impuls or pink noise.
turn into a control signal (average peak or RMS)
Normalize the sound (direct signal = 0 dB)
reverse the control-signal.
search for the different levels with sref.
calculate the RT60 with the gained points (which are higher than the ambient level, i.e. the value for T0)
These calculations should of course be done for all octaves (1/3 octaves) and the labels could be given out for each reverb time / band. .
However, this will most likely look somewhat messy, I fear.
By the way, instead of reversing the sound, it is also possible to work with the reciprikel:
1 (0 dB) → 1; 0.001 (-60 dB) → 1000.0.
