LAFmax from Nyquist plugin to GNU/Octave

[enneric]

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.

Code: Select all

    ;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, [email protected]
    ; 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.

Code: Select all

(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)

[steve]

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

Code: Select all

(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:

Code: Select all

(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/ ... l#index733
snd-avg: http://www.cs.cmu.edu/~rbd/doc/nyquist/ ... l#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.

[enneric]

Thanks for answering steve,

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

Code: Select all

;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, [email protected]
; 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*

[Robert J. H.]

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:

Code: Select all

; 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:

Code: Select all

x=(-1.0:0.01:1.0)

Let's do the biquad section now for Nyquist:

Code: Select all

(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:

Code: Select all

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

Code: Select all

(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

Code: Select all

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

[enneric]

Thank you Rober for your answer,

Yes I have already written some code for this:

Code: Select all

# 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.

Code: Select all

;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, [email protected]
; 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*

[Robert J. H.]

Thank you Rober for your answer,

Yes I have already written some code for this:

Code: Select all

# 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.

Code: Select all

;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, [email protected]
; 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