LAFmax from Nyquist plugin to GNU/Octave

dB LAFmax from Nyquist plugin to GNU/Octave

Hi All,

I’m trying to compute with GNU/Octave, the maximum level with time FAST weighting (LAFmax) from a signal.

On the forum, I found a plugin from mr.ogren that can do it.

    ;nyquist plug-in
    ;version 1
    ;type analyze
    ;name "Equivalent and maximum dB(A)..."
    ;action "Calc. A-weighted equivalent level (LAeq) and maximum level with time weighting FAST (LAFmax)..."

    ; Mikael Ogren, mr.ogren@gmail.com
    ; 2007-01-12
    ; Licensed under GPL, no warranty, use at your own risk...
     
    ; Calibration so that a 1000 Hz tone with amplitude 1.0 gives 94 dB
    (setq calibration (+ 94 28.2))

    ; A-weighting by Edgar (thanks!)
    (setq sa (lp (lp (hp (hp (hp (hp s 20.6) 20.6) 107.7) 737.9) 12200) 12200) )

    ; Exponential time-weighting filter FAST (125 ms)
    ; snd-avg is used to downsample to 100 Hz (by averaging over 441 samples)
    ; This only works for 44.1 kHz sampling frequency, perhaps someone can help out here
    ; by making a more general approach that works for all sampl. frq?
    ; The filtering part is OK for all frequencies, but the "441" constant is not.
    ; -
    ; The constant 0.000001 is to avoid clipping at filtered squared pressure > 1.0
    (setq saf2
    (mult 0.000001
    (snd-avg (snd-biquad (mult sa sa ) 1 0 0 (exp (/ 1 (mult (snd-srate sa) -0.125))) 0 0 0) 441 441 OP-AVERAGE)
    )
    )

    ; Length of the downsampled pressure squared signal
    (setq mlength
    (snd-length saf2 99999999999)
    )

    ; Calc. the equivalent level
    (setq leq
    (+ calibration (* 0.5 (linear-to-db (snd-maxsamp (snd-avg saf2 mlength mlength OP-AVERAGE) ))))
    );

    ; Calc. the maximum level
    (setq lmax
    (+ calibration (* 0.5 (linear-to-db (snd-maxsamp saf2))))
    );

    ; Set the output format to 3 digits (example: 53.3 dB)
    (setq *float-format* "%#3.3g");

    ; Output result as a label track (or append into existing label track)
    (setq u (format NIL "LAeq= ~A  LAFmax= ~A" leq lmax))
    (list (list 0.0 u))

I think that what I need is that line, but I can’t figure out how to convert it in GNU/Octave.

(snd-avg (snd-biquad (mult sa sa ) 1 0 0 (exp (/ 1 (mult (snd-srate sa) -0.125))) 0 0 0) 441 441 OP-AVERAGE)

I don’t use Octave, but perhaps I can help with the Nyquist part.

(snd-avg (snd-biquad (mult sa sa ) 1 0 0 (exp (/ 1 (mult (snd-srate sa) -0.125))) 0 0 0) 441 441 OP-AVERAGE)

There are lots of nested functions here. Let’s separate it out a bit:

(let ((sig (mult sa sa))
      (b0 1)
      (b1 0)
      (b2 0)
      (a1 (exp (/ 1 (mult (snd-srate sa) -0.125))))
      (a2 0)
      (z1init 0)
      (z2init 0))
  (setf sig (snd-biquad sig b0 b1 b2 a1 a2 z1init z2init))
  (snd-avg sig 441 441 OP-AVERAGE))

Here we use a “let” to assign values to local variables: sig, b0, b1, b2, a1, a2, z1init and z2init
We then set the value of the sound “sig” to the result of the biquad evaluation.
Finally we take the average of each 441 samples.

Relevant parts of the Nyquist manual:
snd-biquad: http://www.cs.cmu.edu/~rbd/doc/nyquist/part8.html#index733
snd-avg: http://www.cs.cmu.edu/~rbd/doc/nyquist/part8.html#index663

One problem I note with this code immediately is that “a1” uses the actual sample rate of the sound “sa”, but “snd-avg” seems to assume that the sample rate is 44100 (assuming that the intention is to average over a 0.01 second window). You may want to change that, or restrict your code to a specific sample rate.

Thanks for answering steve,

I found newer version of plugin. It seems to take the actual the sample rate instead of 44100.

;nyquist plug-in
;version 1
;type analyze
;name "Equivalent and maximum dB(A)..."
;action "Calc. A-weighted equivalent level (LAeq) and maximum level with time weighting FAST (LAFmax)..."

; Mikael Ogren, mr.ogren@...
; 2007-01-12
; Licensed under GPL, no warranty, use at your own risk...
;
; Set the output format to 3 digits (example: 53.3 dB)
(setq *float-format* "%#3.3g");
; A-Weighting Filter
(defun cascade (s freqs &aux (f (car freqs)))
  (cond
    ((null f) s)
    ((minusp f) (cascade (lp s (- f)) (cdr freqs)))
    (t (cascade (hp s f) (cdr freqs)))))
;;
;; Main
(let* (
   ; Calibration so that a 1000 Hz tone with amplitude 1.0
   ; gives 94 dB
   (calibration (+ 94 28.2))
   ; A-weighting frequencies by Edgar (thanks!)
   (freqs '(20.6 20.6 107.7 737.9 -12200 -12200))
   (sa (cascade s freqs))
   (step (round (/ *sound-srate* 100.0)))
   (saf2 (snd-avg (biquad (prod sa sa)
                          1e-6   0 0 1   (exp (/ -8 (snd-srate sa))) 0)
                  step step OP-AVERAGE))
   ; Length of the downsampled pressure squared signal
   (mlength(snd-length saf2 ny:all))
   ; Calc. the equivalent level
   (leq (+ calibration (* 0.5 (linear-to-db (snd-maxsamp (snd-avg saf2 mlength mlength OP-AVERAGE))))))
   ; Calc. the maximum level
   (lmax (+ calibration (* 0.5 (linear-to-db (snd-maxsamp saf2)))));
   ; Output result as a label track (or append into existing label track)
   (u (format nil "LAeq= ~A  LAFmax= ~A" leq lmax)))
 u; Sends output to message box if next is commented out
; (list (list 0.0 u))
); end let*

Hi enneric
Do you have already written some code for this?
It seems to me that the a-weighting is a little tricky.
You should design the appropriate 1-order butterworth filter prototypes for “LP” and “HP”.

I’m afraid that the main line you’re interested in can only be translated step by step and controlling the results numerically.

Let’s create a test signal “x” in both languages:

; Nyquist
(setf x (snd-pwl 0 1 '(0 -1.0 200 1.0 201)))
(snd-display x)

Paste the code into the Nyquist prompt and press Debug.

The signal has 201 values from -1.00, -0.99 to 1.00.

In Octave, the same would be:

x=(-1.0:0.01:1.0)

Let’s do the biquad section now for Nyquist:

(setq *float-format* "%#3.5g");
(setf x (snd-pwl 0 1 '(0 -1.0 200 1.0 201)))
(setf x (biquad x 1e-6 0 0 1 (exp (/ -8 (snd-srate x))) 0))
(snd-display x)

Octave uses the Matlab kind of implementation.
This means that A1 and A2 must change sign.

We have x already defined, thus I give here just the filter stage:

x=filter([1e-06, 0, 0], [1, -exp(-8/1), 0], x)

The “/1” should be replaced by “/Fs”, where “Fs” is the predifined sample rate you want to work with eventually.


Let’s skip some steps and produce the output after the averaging in 1 second steps (originally 0.01) at a sample rate of 100 Hz.
The result from this Nyquist code

(setq *float-format* "%#3.5g");
(setf Fs 100)
(setf x (snd-pwl 0 Fs '(0 -1.0 200 1.0 201)))
(setf x (biquad (prod x x) 1e-6 0 0 1 (exp (/ -8 (snd-srate x))) 0))
(snd-display (snd-avg x 100 100 op-average))

should be the same as for the Octave script

Fs=100;
x=(-1.0:0.01:1.0);
x=filter([1e-06, 0, 0], [1, -exp(-8/Fs), 0], x.*x);
newlen=round(length(x)/100)*100
x=postpad(x, newlen);
rtn=sum(reshape(x, Fs, newlen/Fs))/Fs

I hope this helps although I’m no Octave expert.

Robert

Thank you Rober for your answer,

Yes I have already written some code for this:

# Definition of analog A-weighting filter according to IEC/CD 1672.
f1 = 20.598997; 
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
A1000 = 1.9997;

NUMs = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];
DENs = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]); 
DENs = conv(conv(DENs,[1 2*pi*f3]),[1 2*pi*f2]);

# Sampling frequency
Fs = 48000;

# Use the bilinear transformation to get the digital filter. 
[B, A] = bilinear(NUMs,DENs,1/Fs); 

# Load and calibrate the signal
# signal.flac has been convert with Audacity, 
# the format is MAT5 (GNU Octave 2.1 / Matlab 5.0)
load "signal.mat"
signal = wavedata * 206.62;

# Apply A-weighting filter
A_signal = filter(B, A, signal);

# A-weighted equivalent level (LAeq) (dB ref. 2e-5 Pa)
LAeq = 20*log10( sqrt(mean(detrend(A_signal,0).^2)) /2e-5 )

# A-weighted equivalent maximum level with time weighting FAST (LAFmax)
x=filter([1e-6, 0, 0], [1, -exp(-8/Fs), 0], A_signal.**2);

The attached file (signal.flac) can be loaded in GNU/Octave once converted with Audacity in MAT5 format (64 bit float).

I have also calibrated the Audacity plugin, so we are comparing like to like.

;nyquist plug-in
;version 1
;type analyze
;name "Equivalent and maximum dB(A)..."
;action "Calc. A-weighted equivalent level (LAeq) and maximum level with time weighting FAST (LAFmax)..."

; Mikael Ogren, mr.ogren@...
; 2007-01-12
; Licensed under GPL, no warranty, use at your own risk...
;
; Set the output format to 3 digits (example: 53.3 dB)
(setq *float-format* "%#3.3g");
; A-Weighting Filter
(defun cascade (s freqs &aux (f (car freqs)))
  (cond
    ((null f) s)
    ((minusp f) (cascade (lp s (- f)) (cdr freqs)))
    (t (cascade (hp s f) (cdr freqs)))))
;;
;; Main
(let* (
   ; Calibration so that a 1000 Hz tone with amplitude 1.0
   ; gives 94 dB
   (calibration (+ 94 71))
   ; A-weighting frequencies by Edgar (thanks!)
   (freqs '(20.6 20.6 107.7 737.9 -12200 -12200))
   (sa (cascade s freqs))
   (step (round (/ *sound-srate* 100.0)))
   (saf2 (snd-avg (biquad (prod sa sa)
                          1e-6   0 0 1   (exp (/ -8 (snd-srate sa))) 0)
                  step step OP-AVERAGE))
   ; Length of the downsampled pressure squared signal
   (mlength(snd-length saf2 ny:all))
   ; Calc. the equivalent level
   (leq (+ calibration (* 0.5 (linear-to-db (snd-maxsamp (snd-avg saf2 mlength mlength OP-AVERAGE))))))
   ; Calc. the maximum level
   (lmax (+ calibration (* 0.5 (linear-to-db (snd-maxsamp saf2)))));
   ; Output result as a label track (or append into existing label track)
   (u (format nil "LAeq= ~A  LAFmax= ~A" leq lmax)))
 u; Sends output to message box if next is commented out
; (list (list 0.0 u))
); end let*

I see, you’re almost done.
Just three lines more, I guess.

Good luck
Robert