<h2>TL;DR</h2>
<p>I am not aware of any way to predict (or know analytically) exactly what you lose by windowing your signal (which is very signal dependent). Thus, there is no constant factor by which you should multiply in order to compensate for the energy lost by the windowing function.</p>
<hr />
<h2>Energy/Power calculation</h2>
<p>The energy of a signal can be calculated either in the time-domain or the frequency domain (as you have already mentioned Parseval's theorem). Since the Fourier Transform is a linear operator, you retain all information when you move from one domain to the other.</p>
<p>The energy of a signal is given by the squared sum of its time-domain samples. There is a mathematical &quot;discrepancy&quot; when it comes to the frequency-domain representation though. When you perform a forward and then an inverse transform you have to divide by the length of your signal vector (i.e. the number of samples you use). Nevertheless, from a mathematical viewpoint, there is nothing wrong if you perform this division at either of the transformations - Forward and Inverse - (dividing with the number of samples/frequency bins) or even at both (dividing by the square root of the number of samples in each transformation).</p>
<p>Many software perform the division at the inverse transformation. This means that you may have to perform the division yourself if you want to get accurate amplitude representation. But, please note that this is implementation dependent. Additionally, if you perform the division to get accurate magnitude values, you may have to multiply by the same number before performing an inverse transform (since the software may perform the division again).</p>
<p>The power of the signal is given by the formula</p>
<pre><code>e = (1/N) * sum(x .* x)
</code></pre>
<p>where the symbol <code>.*</code> denotes element-wise multiplication (or else the Hadamard product). This formula holds both for the frequency and the time domains. Now, of course this is without the application of a windowing function to your signal (identical to using the rectangular window). When you apply a window you lose information and keep in mind that even without taking into account the numerical precision, this process may be irreversible, when the window goes to zero at the edges.</p>
<h2>Inspiration from COLA</h2>
<p>Now, looking further into the <a href="https://www.dsprelated.com/freebooks/sasp/Constant_Overlap_Add_COLA_Cases.html" rel="nofollow noreferrer"><strong>Constant Overlap Add Constraint</strong> (<strong>COLA</strong>)</a> we see that for a <em>Hann</em> window with overlap equal to 50% you can have perfect reconstruction, which means that you don't lose any information of the signal. Similar COLA constraints can be found in the frequency domain (see <a href="https://www.dsprelated.com/freebooks/sasp/Frequency_Domain_COLA_Constraints.html" rel="nofollow noreferrer">this</a>). Now, one could assume that since a 50% overlap provides all the information, a Hann window halves the frequency of the signal. This is <strong>not</strong> true though. You can have simple experiments with random signals (Gaussian or uniformly distributed signals are adequate) and calculate the energy in the time domain and the windowed signal's energy and you'll see that there is not a constant relation but the results are very signal-dependent. One remark here is that COLA constraint (as well as other similar constraints) refer to overlapping segments of a signal, so in this way what you lose from one frame will be present in the next one and this is why you retain all information in such a (usually STFT) scheme.</p>
<h2>Example</h2>
<p>A simple intuitive example to realise that you may lose information is to think of an impulsive sound (like clicks I guess but I am not very familiar with echo-location signals or other bioacoustics techniques and terms) which is located at the beginning of the signal. By windowing you'll attenuate a big portion of the energy. On the contrary, if you centre the impulse you'll lose but a small fraction (depending on the type and length of the window of course).</p>

CC BY-SA 4.0
2022-10-09T00:05:55.520Z
2022-10-09T00:05:55.520Z
ZaellixA
<p>As far I understand this question, there is a lot of confusion in terminology. So, I tried to summarize my understanding.
All methods below are based on the magnitude squared output of the FFT: p=|FFT(x)|^2</p>
<p>Depending on the quantity you wanted to estimate you apply different scaling</p>
<ul>
<li><p>power spectral density (PSD) of noise</p>
<blockquote>
<p>Scale FFT power estimate with $1/(f_s \sum w^2[n])$ to obtain power/Hz</p>
</blockquote>
</li>
<li><p>RMS estimates of tonals</p>
<blockquote>
<p>Scale FFT peak power estimate with $1/(\sum w[n])^2$ to obtain signal rms
estimate</p>
</blockquote>
</li>
<li><p>peak amplitude estimates of transients (echolocation clicks)</p>
<blockquote>
<p>Scale FFT peak power estimate with $1/N$ to obtain signal peak amplitude</p>
</blockquote>
</li>
</ul>
<p>If someone wanted to play with my test program (python jupyter notebook) <a href="https://wmxz-eu.github.io/power_spectra.html" rel="nofollow noreferrer">here is the link</a>. Caveat: work in progress</p>
<p>Edit: to address the comment &quot;how to correctly scale a power spectrum (not peak amp estimate) for an echolocation click when a window function&quot;, we have to make clear scale to read what? I tried to demonstrate (with the jupyter notebook (link above)) that both methods 'density' and 'spectrum' estimates (in Python language) do not make sense for echolocation clicks that are much shorted than analysis window.</p>
<p>Spectral analysis of echolocation clicks is made to see the shape of the spectral distribution. While most publication recognize this and present the spectral distribution normalized to peak level, a possible (and IMHO the only correct) scaling is to normalize so that spectral peak level corresponds to timeseries peak level.</p>
<p>Scaling to power/Hz or RMS gives the impression that echolocation clicks are continuous signals where noise parameters (power/Hz) or tonal parameters (RMS) make sense.</p>

CC BY-SA 4.0
2022-10-11T22:35:16.920Z
This is helpful but doesn't directly answer my question for how to correctly scale a power spectrum (not peak amp estimate) for an echolocation click when a window function (e.g., Hann) is used. Could you add to your answer to include this?
Marie
CC BY-SA 4.0
2022-10-12T16:05:29.717Z
Thank you editing your response. Your developing jupyter notebook is and will be a great resource.
Marie
CC BY-SA 4.0
2022-10-11T20:07:22.400Z
2022-10-12T05:15:06.077Z
2022-10-12T05:15:06.077Z
WMXZ
<p>In my ongoing research, I found an answer by <a href="https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.12330" rel="nofollow noreferrer">Merchant et al. (2015)</a>. There, the authors summarize acoustic measurements in detail. This paper is also referenced in <a href="https://link.springer.com/book/10.1007/978-3-030-97540-1" rel="nofollow noreferrer">this textbook</a>. In the supplemental material, there is a PDF tutorial (<a href="https://besjournals.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1111%2F2041-210X.12330&amp;file=mee312330-sup-0001-AppendixS1.pdf" rel="nofollow noreferrer">link here</a>) that explains all calculations conducted in PAMGuide and includes complimentary R code (PAMGuide.zip). A description of how PSDs are calculated begins on page 29 of the tutorial.</p>
<p>If you read through that document and look at the code for PSD, there is one additional scaling step that they include at the end which accounts for the frequency resolution (i.e., how many Hz per bin determined by the sampling rate and FFT size). Below I revised my example code to match the PSD steps used by Merchant et al. (2015). The authors distinguish between 'power spectrum' and 'power spectral densities', where the former does not have the extra scaling step to standardize by 1 Hz bins.</p>
<pre><code># signal x
x &lt;- runif(512)
N &lt;- length(x) # FFT length

# Multiply by Hann window function
window &lt;- signal::hanning(N)
sig_wind &lt;- x * window
# Scale result to account for lost amplitude due to the window function
# alpha=0.5 in the PAMGuide tutorial which is equal to mean(window) for Hann
sig_wind &lt;- sig_wind / mean(window)

# perform FFT and convert to intensity (power)
# divide by N here
power &lt;- (Mod(fft(sig_wind)) / N)^2

# the power spectrum ('power') is symmetrical about the Nyquist freq, 
# so frequencies above this are discarded
# remaining frequencies are doubled to account for total energy 
# transformed (Parseval’s theorem)
power_spec &lt;- 2 * power[0:round(N/2)]

# convert to dB and standardize to 1 Hz bin widths
# this is equation 11 in the tutorial and in the R code
# (resulting units of dB re 1 uPa^2/Hz)
# from the tutorial: 
# “B is the noise power bandwidth of the window function, 
# which corrects for the energy added through spectral leakage”
B &lt;- 1.5 # for Hann window
# freq resolution = sampling rate divided by N (i.e., how many Hz for each bin)
delt_freq &lt;- samp_freq / N
scale &lt;- 1 / B * delt_freq
PSD_dB &lt;- 10*log10(power_spec/p_ref^2 * scale)# p_ref^2 = 1
</code></pre>
<p>However, this PSD calculation differs from others (e.g., see <a href="https://www.oceansinitiative.org/wp-content/uploads/2012/07/PocketBook-3rd-ed.pdf" rel="nofollow noreferrer">here</a>). Perhaps there are simply many approaches to this problem of PSD estimations and scaling (evidenced by <a href="https://dsp.stackexchange.com/questions/2096/why-so-many-methods-of-computing-psd">this thread</a>).</p>

CC BY-SA 4.0
2022-10-11T17:56:14.787Z
2022-10-11T17:56:14.787Z
Marie
<p>I am plotting power spectra of echolocation clicks in R and I want to calculate the absolute magnitudes of the power spectra (not normalized dB values relative to the max). My understanding is that there can be a scaling factor applied in two possible steps when calculating the PSD of a click:</p>
<ol>
<li>when a window function is applied (e.g., hanning window)</li>
<li>to the FFT output</li>
</ol>
<p>For #1: A similar question regarding scaling prior to the fft() function was posted recently <a href="https://bioacoustics.stackexchange.com/questions/1179/why-divide-by-average-of-window-when-calculating-spectrum-of-signal">here</a>. In this case, a hanning function was applied to a signal and the result could be scaled by the mean(window), mean(window^2), or sum(window^2).</p>
<p>For #2: The R documentation for the fft() function states: &quot;the value computed and returned by fft is the <em>unnormalized</em> univariate discrete Fourier transform of the sequence of values in z.&quot; There are a variety of explanations for scaling/normalizing the output of FFT functions (e.g., <a href="https://electronics.stackexchange.com/questions/25900/scaling-fft-output-by-number-of-points-in-fft">here</a> and <a href="https://dsp.stackexchange.com/questions/66058/am-i-supposed-to-normalize-fft-in-python">here</a>) where some multiply the output by 1/N or 1/sqrt(N), where N is the FFT size. To test whether Parseval's theorem of equal energy is satisfied, I compared the squared sum of the waveform values and the squared sum of the spectrum. For these to be equal, I had to normalize the fft() by 1/sqrt(N).</p>
<p>Building on the code referenced above (found <a href="https://bioacoustics.stackexchange.com/questions/1179/why-divide-by-average-of-window-when-calculating-spectrum-of-signal">here</a>):</p>
<pre><code># signal x
x &lt;- runif(512)
window &lt;- signal::hanning(length(x))
sig_wind &lt;- x * window
# why this step?? or should it be mean(window^2) ? or sum(window^2)??
sig_wind &lt;- sig_wind / mean(window)

# compute fft without scaling
spec &lt;- 20*log10(abs(fft(sig_wind)))
# or compute fft and normalize it by the FFT size or /sqrt(length(sig_wind))
spec &lt;- 20*log10(abs(fft(sig_wind))/length(sig_wind))
# ^ which is correct?
</code></pre>
<p>What is the correct approach for scaling/normalizing the input (if a window function is applied) and output of an FFT so that the energy calculated for the power spectrum accurately reflects that in the waveform?</p>

CC BY-SA 4.0
2022-10-06T13:04:20.213Z
Thanks for posting your question on BSE! Can you spell out ‘PSD’ for clarity?
Shannon
CC BY-SA 4.0
2022-10-06T16:58:02.190Z
Yes! Just edited.
Marie
CC BY-SA 4.0
2022-10-10T15:51:39.843Z
You cannot estimate the PSD for echolocation clicks, as PSD makes only sense for random processes and PSD does never reflect the waveform. I put this as comment, as any correct answer will not be helpful for the question asked.
WMXZ
CC BY-SA 4.0
2022-10-10T17:57:15.370Z
@WMXZ In bioacoustics, power spectra are frequently computed for echolocation clicks. [Madsen and Wahlberg (2007)](https://www.sciencedirect.com/science/article/pii/S0967063707001124) provide a comprehensive overview and here are a few other papers as examples where this is applied: [here](https://asa.scitation.org/doi/10.1121/1.4817832), [here](https://www.navymarinespeciesmonitoring.us/files/7913/8308/8904/Baumann-Pickering_et_al._2013-Aleutian_Islands_beaked_whale_echolocation_signals.pdf), and [here](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0162069)
Marie
CC BY-SA 4.0
2022-10-10T18:02:58.500Z
@WMXZ Can you please explain further why a PSD doesn't make sense?
Marie
CC BY-SA 4.0
2022-10-10T18:16:02.917Z
power spectra YES, but power spectra density NO (i.e. uPa^2/Hz) for this you need either noise as long as FFT or length of click.  Also, PSD is defined via the autocorrelation function that you cannot obtain from a click. If the term PSD is only a typo, please rewrite question. See also https://dsp.stackexchange.com/questions/4691/what-is-the-difference-between-psd-and-squared-magnitude-of-frequency-spectrum?lq=1
WMXZ
CC BY-SA 4.0
2022-10-11T17:38:31.137Z
@WMXZ It is common to see normalized power spectra reported in bioacoustics papers (implied units of dB re 1 uPa^2) but I also see units of 're 1 uPa^2/Hz' for spectra of echolocation clicks (e.g., [here](https://link.springer.com/article/10.1007/s00300-022-03008-5)). Are you saying those units are incorrectly reported or am I missing something? My understanding is that there are many ways to estimate PSD (see [here](https://dsp.stackexchange.com/questions/2096/why-so-many-methods-of-computing-psd)). thanks for your thoughts.
Marie
<signal-processing><fft><window><power-spectral-density>
How and why do we rescale the input and/or output of an FFT for bioacoustics research?
# How and why do we rescale the input and/or output of an FFT for bioacoustics research?

(Asked by: username_0 on 2022-10-05)

I am plotting power spectra of echolocation clicks in R and I want to calculate the absolute magnitudes of the power spectra (not normalized dB values relative to the max). My understanding is that there can be a scaling factor applied in two possible steps when calculating the PSD of a click:

1. when a window function is applied (e.g., hanning window)
2. to the FFT output

For #1: A similar question regarding scaling prior to the fft() function was posted recently [here](https://bioacoustics.stackexchange.com/questions/1179/why-divide-by-average-of-window-when-calculating-spectrum-of-signal). In this case, a hanning function was applied to a signal and the result could be scaled by the mean(window), mean(window^2), or sum(window^2).

For #2: The R documentation for the fft() function states: "the value computed and returned by fft is the *unnormalized* univariate discrete Fourier transform of the sequence of values in z." There are a variety of explanations for scaling/normalizing the output of FFT functions (e.g., [here](https://electronics.stackexchange.com/questions/25900/scaling-fft-output-by-number-of-points-in-fft) and [here](https://dsp.stackexchange.com/questions/66058/am-i-supposed-to-normalize-fft-in-python)) where some multiply the output by 1/N or 1/sqrt(N), where N is the FFT size. To test whether Parseval's theorem of equal energy is satisfied, I compared the squared sum of the waveform values and the squared sum of the spectrum. For these to be equal, I had to normalize the fft() by 1/sqrt(N).

Building on the code referenced above (found [here](https://bioacoustics.stackexchange.com/questions/1179/why-divide-by-average-of-window-when-calculating-spectrum-of-signal)):

```
# signal x
x <- runif(512)
window <- signal::hanning(length(x))
sig_wind <- x * window
# why this step?? or should it be mean(window^2) ? or sum(window^2)??
sig_wind <- sig_wind / mean(window)

# compute fft without scaling
spec <- 20*log10(abs(fft(sig_wind)))
# or compute fft and normalize it by the FFT size or /sqrt(length(sig_wind))
spec <- 20*log10(abs(fft(sig_wind))/length(sig_wind))
# ^ which is correct?

```

What is the correct approach for scaling/normalizing the input (if a window function is applied) and output of an FFT so that the energy calculated for the power spectrum accurately reflects that in the waveform?

**Tags:** signal-processing, fft, window, power-spectral-density

> **Comments:**
> **username_1**: Thanks for posting your question on BSE! Can you spell out ‘PSD’ for clarity?  
> **username_0**: Yes! Just edited.  
> **username_2**: You cannot estimate the PSD for echolocation clicks, as PSD makes only sense for random processes and PSD does never reflect the waveform. I put this as comment, as any correct answer will not be helpful for the question asked.  
> **username_0**: @username_2 In bioacoustics, power spectra are frequently computed for echolocation clicks. [Madsen and Wahlberg (2007)](https://www.sciencedirect.com/science/article/pii/S0967063707001124) provide a comprehensive overview and here are a few other papers as examples where this is applied: [here](https://asa.scitation.org/doi/10.1121/1.4817832), [here](https://www.navymarinespeciesmonitoring.us/files/7913/8308/8904/Baumann-Pickering_et_al._2013-Aleutian_Islands_beaked_whale_echolocation_signals.pdf), and [here](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0162069)  
> **username_0**: @username_2 Can you please explain further why a PSD doesn't make sense?  
> **username_2**: power spectra YES, but power spectra density NO (i.e. uPa^2/Hz) for this you need either noise as long as FFT or length of click.  Also, PSD is defined via the autocorrelation function that you cannot obtain from a click. If the term PSD is only a typo, please rewrite question. See also https://dsp.stackexchange.com/questions/4691/what-is-the-difference-between-psd-and-squared-magnitude-of-frequency-spectrum?lq=1  
> **username_0**: @username_2 It is common to see normalized power spectra reported in bioacoustics papers (implied units of dB re 1 uPa^2) but I also see units of 're 1 uPa^2/Hz' for spectra of echolocation clicks (e.g., [here](https://link.springer.com/article/10.1007/s00300-022-03008-5)). Are you saying those units are incorrectly reported or am I missing something? My understanding is that there are many ways to estimate PSD (see [here](https://dsp.stackexchange.com/questions/2096/why-so-many-methods-of-computing-psd)). thanks for your thoughts.  

---

## 3 answers

### Answer by username_3 on 2022-10-09. Score: 4

## TL;DR

I am not aware of any way to predict (or know analytically) exactly what you lose by windowing your signal (which is very signal dependent). Thus, there is no constant factor by which you should multiply in order to compensate for the energy lost by the windowing function.

---

## Energy/Power calculation

The energy of a signal can be calculated either in the time-domain or the frequency domain (as you have already mentioned Parseval's theorem). Since the Fourier Transform is a linear operator, you retain all information when you move from one domain to the other.

The energy of a signal is given by the squared sum of its time-domain samples. There is a mathematical "discrepancy" when it comes to the frequency-domain representation though. When you perform a forward and then an inverse transform you have to divide by the length of your signal vector (i.e. the number of samples you use). Nevertheless, from a mathematical viewpoint, there is nothing wrong if you perform this division at either of the transformations - Forward and Inverse - (dividing with the number of samples/frequency bins) or even at both (dividing by the square root of the number of samples in each transformation).

Many software perform the division at the inverse transformation. This means that you may have to perform the division yourself if you want to get accurate amplitude representation. But, please note that this is implementation dependent. Additionally, if you perform the division to get accurate magnitude values, you may have to multiply by the same number before performing an inverse transform (since the software may perform the division again).

The power of the signal is given by the formula

```
e = (1/N) * sum(x .* x)

```

where the symbol `.*` denotes element-wise multiplication (or else the Hadamard product). This formula holds both for the frequency and the time domains. Now, of course this is without the application of a windowing function to your signal (identical to using the rectangular window). When you apply a window you lose information and keep in mind that even without taking into account the numerical precision, this process may be irreversible, when the window goes to zero at the edges.

## Inspiration from COLA

Now, looking further into the [**Constant Overlap Add Constraint** (**COLA**)](https://www.dsprelated.com/freebooks/sasp/Constant_Overlap_Add_COLA_Cases.html) we see that for a *Hann* window with overlap equal to 50% you can have perfect reconstruction, which means that you don't lose any information of the signal. Similar COLA constraints can be found in the frequency domain (see [this](https://www.dsprelated.com/freebooks/sasp/Frequency_Domain_COLA_Constraints.html)). Now, one could assume that since a 50% overlap provides all the information, a Hann window halves the frequency of the signal. This is **not** true though. You can have simple experiments with random signals (Gaussian or uniformly distributed signals are adequate) and calculate the energy in the time domain and the windowed signal's energy and you'll see that there is not a constant relation but the results are very signal-dependent. One remark here is that COLA constraint (as well as other similar constraints) refer to overlapping segments of a signal, so in this way what you lose from one frame will be present in the next one and this is why you retain all information in such a (usually STFT) scheme.

## Example

A simple intuitive example to realise that you may lose information is to think of an impulsive sound (like clicks I guess but I am not very familiar with echo-location signals or other bioacoustics techniques and terms) which is located at the beginning of the signal. By windowing you'll attenuate a big portion of the energy. On the contrary, if you centre the impulse you'll lose but a small fraction (depending on the type and length of the window of course).


---

### Answer by username_2 on 2022-10-11. Score: 3

As far I understand this question, there is a lot of confusion in terminology. So, I tried to summarize my understanding.
All methods below are based on the magnitude squared output of the FFT: p=|FFT(x)|^2

Depending on the quantity you wanted to estimate you apply different scaling

* power spectral density (PSD) of noise

  > Scale FFT power estimate with $1/(f\_s \sum w^2[n])$ to obtain power/Hz
* RMS estimates of tonals

  > Scale FFT peak power estimate with $1/(\sum w[n])^2$ to obtain signal rms
  > estimate
* peak amplitude estimates of transients (echolocation clicks)

  > Scale FFT peak power estimate with $1/N$ to obtain signal peak amplitude

If someone wanted to play with my test program (python jupyter notebook) [here is the link](https://username_2-eu.github.io/power_spectra.html). Caveat: work in progress

Edit: to address the comment "how to correctly scale a power spectrum (not peak amp estimate) for an echolocation click when a window function", we have to make clear scale to read what? I tried to demonstrate (with the jupyter notebook (link above)) that both methods 'density' and 'spectrum' estimates (in Python language) do not make sense for echolocation clicks that are much shorted than analysis window.

Spectral analysis of echolocation clicks is made to see the shape of the spectral distribution. While most publication recognize this and present the spectral distribution normalized to peak level, a possible (and IMHO the only correct) scaling is to normalize so that spectral peak level corresponds to timeseries peak level.

Scaling to power/Hz or RMS gives the impression that echolocation clicks are continuous signals where noise parameters (power/Hz) or tonal parameters (RMS) make sense.

> **Comments:**
> **username_0**: This is helpful but doesn't directly answer my question for how to correctly scale a power spectrum (not peak amp estimate) for an echolocation click when a window function (e.g., Hann) is used. Could you add to your answer to include this?  
> **username_0**: Thank you editing your response. Your developing jupyter notebook is and will be a great resource.  

---

### Answer by username_0 on 2022-10-11. Score: 1

In my ongoing research, I found an answer by [Merchant et al. (2015)](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.12330). There, the authors summarize acoustic measurements in detail. This paper is also referenced in [this textbook](https://link.springer.com/book/10.1007/978-3-030-97540-1). In the supplemental material, there is a PDF tutorial ([link here](https://besjournals.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1111%2F2041-210X.12330&file=mee312330-sup-0001-AppendixS1.pdf)) that explains all calculations conducted in PAMGuide and includes complimentary R code (PAMGuide.zip). A description of how PSDs are calculated begins on page 29 of the tutorial.

If you read through that document and look at the code for PSD, there is one additional scaling step that they include at the end which accounts for the frequency resolution (i.e., how many Hz per bin determined by the sampling rate and FFT size). Below I revised my example code to match the PSD steps used by Merchant et al. (2015). The authors distinguish between 'power spectrum' and 'power spectral densities', where the former does not have the extra scaling step to standardize by 1 Hz bins.

```
# signal x
x <- runif(512)
N <- length(x) # FFT length

# Multiply by Hann window function
window <- signal::hanning(N)
sig_wind <- x * window
# Scale result to account for lost amplitude due to the window function
# alpha=0.5 in the PAMGuide tutorial which is equal to mean(window) for Hann
sig_wind <- sig_wind / mean(window)

# perform FFT and convert to intensity (power)
# divide by N here
power <- (Mod(fft(sig_wind)) / N)^2

# the power spectrum ('power') is symmetrical about the Nyquist freq, 
# so frequencies above this are discarded
# remaining frequencies are doubled to account for total energy 
# transformed (Parseval’s theorem)
power_spec <- 2 * power[0:round(N/2)]

# convert to dB and standardize to 1 Hz bin widths
# this is equation 11 in the tutorial and in the R code
# (resulting units of dB re 1 uPa^2/Hz)
# from the tutorial: 
# “B is the noise power bandwidth of the window function, 
# which corrects for the energy added through spectral leakage”
B <- 1.5 # for Hann window
# freq resolution = sampling rate divided by N (i.e., how many Hz for each bin)
delt_freq <- samp_freq / N
scale <- 1 / B * delt_freq
PSD_dB <- 10*log10(power_spec/p_ref^2 * scale)# p_ref^2 = 1

```

However, this PSD calculation differs from others (e.g., see [here](https://www.oceansinitiative.org/wp-content/uploads/2012/07/PocketBook-3rd-ed.pdf)). Perhaps there are simply many approaches to this problem of PSD estimations and scaling (evidenced by [this thread](https://dsp.stackexchange.com/questions/2096/why-so-many-methods-of-computing-psd)).


---

