Answers.Body: <p>In order to demonstrate the relation between noise spectral density, noise bandwidth and rms estimation I offer a slightly modified example code</p>
  <pre><code>library(tuneR)
  library(soundgen)
  library(seewave)
  #
  noise_power=1e-3                  # spectral noise power V^2/Hz
  samp_frq=24000                    # sampling frequency
  noise_BW = samp_frq/2             # noise bandwidth
  noise_var=noise_power*noise_BW    # noise variance V^2/Hz * Hz 
   
  #
  # noise time series
  noise = rnorm(5*samp_frq,sd=sqrt(noise_var))
  #
  # convert to wav object
  w = Wave(left = noise, samp.rate = samp_frq, bit=32)
  
  # custom spectral averaging
  f &lt;- function(x){  return (mean(abs(x)^2)) }
  #
  # processing
  window_length_in_ms=100
  window_length_in_samples=(w@samp.rate * window_length_in_ms) / 1000
  #
  rms = getRMS(w, windowLength = window_length_in_ms, 
               normalize = F,  overlap = 50)
  power_spec = meanspec(w, wn=&quot;rectangle&quot;, wl = window_length_in_samples, 
                        FUN=f, 
                        norm=F, ovlp = 50, plot = F)
  #
  # print results
  print(10*log10(noise_power))
  print(20*log10(sd(noise))               -10*log10(noise_BW))
  print(20*log10(rms<span class="math-container">$summary$</span>ampl_mean)   -10*log10(noise_BW))
  print(10*log10(sum(power_spec[,&quot;y&quot;])/2) -10*log10(noise_BW))
  </code></pre>
  <p>the results of which are</p>
  <pre><code>[1] -30
  [1] -30.01798
  [1] -30.01745
  [1] -30.01693
  </code></pre>
  <p>The users are invited to play with sampling frequency and window length to note that the estimated (rms-based) noise power levels remain constant at -30dB or 1e-3</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-12-12T16:53:37.457Z
Answers.LastActivityDate: 2023-12-14T06:21:39.423Z
Answers.LastEditDate: 2023-12-14T06:21:39.423Z
Answers.UserDisplayName: WMXZ
Answers.Body: <p>The differences are due to different scalings.</p>
  <p>the following modifications give all the same answer (or similar due to windowing)</p>
  <pre><code># now everything relative to amplitude unit 1
  f &lt;- function(x)
  {
    return(mean(abs(x)^2))
  }
  
  std = sqrt(mean(w@left^2))
  
  rms = getRMS(w, windowLength = window_length_in_ms, normalize = F,  overlap = 50)
  power_spec = meanspec(w, wn=&quot;rectangle&quot;, wl = window_length_in_samples, norm=F, FUN=f, ovlp = 50, plot = T)
  
  print(20*log10(1000))
  print(20*log10(std))
  print(20*log10(rms<span class="math-container">$summary$</span>ampl_mean))
  print(10*log10(sum(power_spec[,&quot;y&quot;])/2))
  </code></pre>
  <p>generating as output</p>
  <pre><code>[1] 60
  [1] 59.98897
  [1] 59.97903
  [1] 59.99018
  </code></pre>
  <p>The obtained level is 60 dB (or similar due to random samples) relative to 1 amplitude unit as it should be.</p>
  <p>The are multiple differences:</p>
  <ol>
  <li>the meanSpec routine seems not to correct for the windowing shape (using different windows gives different results) So, rectangle window should be used to compare</li>
  <li>the SeeWave meanSpec function uses as default average the mean of the modulus of the spectrum. This is physically not correct. The power (or modulus square) should be used for averaging. Therefore I added a custom function.</li>
  <li>To convert the power spectrum to dB, 10*log10(.) must be used</li>
  <li>to convert power spectrum to broadband rms values, the sum over all
  frequencies should be used and not the mean. This can easily be verified by changing the window length.The factor 2 is still unclear and most likely due to the stdft function.</li>
  </ol>
  <p>I would strongly recommend to use scaling to max value (in data or spectrum) ONLY for visualization purpose, and otherwise scale to 1 (or acoustic reference 1uPa or 20uPa if sensitivity is known).</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-12-12T12:03:02.467Z
Answers.LastActivityDate: 2023-12-12T12:03:02.467Z
Answers.UserDisplayName: WMXZ
Body: <p>I want to calculate the sound pressure level (in dB) of a wave and then display the mean power spectrum.  The issue I face is that the SPL value I get from getRMS() is very different from the mean amplitude values from meanspec().</p>
  <pre><code>library(tuneR)
  library(soundgen)
  library(seewave)
  
  #make random noise
  noise = rnorm(24000*5, sd = 1000)
  w = Wave(left = noise, samp.rate = 24000)
  
  #set window lengths (ms for getRMS, samples for meanspec)
  window_length_in_ms = 50
  window_length_in_samples = w@samp.rate * (window_length_in_ms / 1000)
  
  #calculate SPL
  rms = getRMS(w, windowLength = window_length_in_ms, overlap = 50)
  SPLdB = ((20 * log10(rms<span class="math-container">$summary$</span>ampl_mean)))
  SPLdB  #this is the same RMS value I get from Audacity
  
  #graph power spectrum
  power_spec = meanspec(w, wl = window_length_in_samples, dB = &quot;max0&quot;, ovlp = 50, plot = T)
  mean(power_spec[,&quot;y&quot;]). #this value is VERY different from my SPL
  </code></pre>
  <p>Obviously, I am missing something fundamental.  I appreciate any feedback.</p>
  
Tags: <r><rms><seewave><soundgen>
Title: Can't reconcile amplitude from getRMS (soundgen) with output of meanspec (seewave) in R
ThreadText: # Can't reconcile amplitude from getRMS (soundgen) with output of meanspec (seewave) in R
  
  (Asked by: username_0 on 2023-12-11)
  
  I want to calculate the sound pressure level (in dB) of a wave and then display the mean power spectrum.  The issue I face is that the SPL value I get from getRMS() is very different from the mean amplitude values from meanspec().
  library(tuneR)
  library(soundgen)
  library(seewave)
  
  #make random noise
  noise = rnorm(24000*5, sd = 1000)
  w = Wave(left = noise, samp.rate = 24000)
  
  #set window lengths (ms for getRMS, samples for meanspec)
  window_length_in_ms = 50
  window_length_in_samples = w@samp.rate * (window_length_in_ms / 1000)
  
  #calculate SPL
  rms = getRMS(w, windowLength = window_length_in_ms, overlap = 50)
  SPLdB = ((20 * log10(rms$summary$ampl_mean)))
  SPLdB  #this is the same RMS value I get from Audacity
  
  #graph power spectrum
  power_spec = meanspec(w, wl = window_length_in_samples, dB = "max0", ovlp = 50, plot = T)
  mean(power_spec[,"y"]). #this value is VERY different from my SPL
  
  Obviously, I am missing something fundamental.  I appreciate any feedback.
  
  **Tags:** r, rms, seewave, soundgen
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2023-12-12. Score: 0 | Accepted answer
  
  In order to demonstrate the relation between noise spectral density, noise bandwidth and rms estimation I offer a slightly modified example code
  library(tuneR)
  library(soundgen)
  library(seewave)
  #
  noise_power=1e-3                  # spectral noise power V^2/Hz
  samp_frq=24000                    # sampling frequency
  noise_BW = samp_frq/2             # noise bandwidth
  noise_var=noise_power*noise_BW    # noise variance V^2/Hz * Hz 
   
  #
  # noise time series
  noise = rnorm(5*samp_frq,sd=sqrt(noise_var))
  #
  # convert to wav object
  w = Wave(left = noise, samp.rate = samp_frq, bit=32)
  
  # custom spectral averaging
  f <- function(x){  return (mean(abs(x)^2)) }
  #
  # processing
  window_length_in_ms=100
  window_length_in_samples=(w@samp.rate * window_length_in_ms) / 1000
  #
  rms = getRMS(w, windowLength = window_length_in_ms, 
               normalize = F,  overlap = 50)
  power_spec = meanspec(w, wn="rectangle", wl = window_length_in_samples, 
                        FUN=f, 
                        norm=F, ovlp = 50, plot = F)
  #
  # print results
  print(10*log10(noise_power))
  print(20*log10(sd(noise))               -10*log10(noise_BW))
  print(20*log10(rms$summary$ampl_mean)   -10*log10(noise_BW))
  print(10*log10(sum(power_spec[,"y"])/2) -10*log10(noise_BW))
  
  the results of which are
  [1] -30
  [1] -30.01798
  [1] -30.01745
  [1] -30.01693
  
  The users are invited to play with sampling frequency and window length to note that the estimated (rms-based) noise power levels remain constant at -30dB or 1e-3
  
  
  ---
  
  ### Answer by username_1 on 2023-12-12. Score: 0
  
  The differences are due to different scalings.
  the following modifications give all the same answer (or similar due to windowing)
  # now everything relative to amplitude unit 1
  f <- function(x)
  {
    return(mean(abs(x)^2))
  }
  
  std = sqrt(mean(w@left^2))
  
  rms = getRMS(w, windowLength = window_length_in_ms, normalize = F,  overlap = 50)
  power_spec = meanspec(w, wn="rectangle", wl = window_length_in_samples, norm=F, FUN=f, ovlp = 50, plot = T)
  
  print(20*log10(1000))
  print(20*log10(std))
  print(20*log10(rms$summary$ampl_mean))
  print(10*log10(sum(power_spec[,"y"])/2))
  
  generating as output
  [1] 60
  [1] 59.98897
  [1] 59.97903
  [1] 59.99018
  
  The obtained level is 60 dB (or similar due to random samples) relative to 1 amplitude unit as it should be.
  The are multiple differences:
  
  the meanSpec routine seems not to correct for the windowing shape (using different windows gives different results) So, rectangle window should be used to compare
  the SeeWave meanSpec function uses as default average the mean of the modulus of the spectrum. This is physically not correct. The power (or modulus square) should be used for averaging. Therefore I added a custom function.
  To convert the power spectrum to dB, 10*log10(.) must be used
  to convert power spectrum to broadband rms values, the sum over all
  frequencies should be used and not the mean. This can easily be verified by changing the window length.The factor 2 is still unclear and most likely due to the stdft function.
  
  I would strongly recommend to use scaling to max value (in data or spectrum) ONLY for visualization purpose, and otherwise scale to 1 (or acoustic reference 1uPa or 20uPa if sensitivity is known).
  
  
  ---
  
  