ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 5 Pull Requests Actions Releases 14 Activity

Restored FFT scaling to power/bin instead of PSD

Browse Source
This commit is contained in:
Casper Jansen 2021-10-22 13:36:23 +02:00
parent ea0edd2817
commit 9c8c546566
2 changed files with 9 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
lasp/c/lasp_ps.c
View File

@ -117,11 +117,15 @@ void PowerSpectra_compute(const PowerSpectra* ps,
TRACE(15,"fft done"); TRACE(15,"fft done");
/* Scale fft such that power is easily computed /* Scale fft such that power is easily comxputed */
const c scale_fac = d_sqrt(2/win_pow)/nfft;
/* POWER SPECTRAL DENSITY? Scale fft such that power is easily computed
* - Multiply power spectral density by 2 except at f=0 and f=fNq * - Multiply power spectral density by 2 except at f=0 and f=fNq
* - Divide by energy of window function = nfft * window_power * - Divide by energy of window function = nfft * window_power
* - .. sqrt(factors) because it is applied to output fft instead of psd */ * - .. sqrt(factors) because it is applied to output fft instead of psd */
const c scale_fac = d_sqrt(2/(nfft*win_pow)); // const c scale_fac = d_sqrt(2/(nfft*win_pow));
cmat_scale(&fft_work,scale_fac); cmat_scale(&fft_work,scale_fac);
TRACE(15,"scale done"); TRACE(15,"scale done");

5
lasp/tools/tools.py
View File

@ -48,6 +48,7 @@ class SmoothingWidth(Enum):
twelve = (12, '1/12th octave smoothing') twelve = (12, '1/12th octave smoothing')
twfo = (24, '1/24th octave smoothing') twfo = (24, '1/24th octave smoothing')
ftei = (48, '1/48th octave smoothing') ftei = (48, '1/48th octave smoothing')
# hundred = (100, '1/100th octave smoothing') # useful for removing 'grass'
@staticmethod @staticmethod
def fillComboBox(cb): def fillComboBox(cb):
@ -136,7 +137,6 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
x0 = 1 if freq[0] == 0 else 0 # Skip first data point if zero frequency x0 = 1 if freq[0] == 0 else 0 # Skip first data point if zero frequency
Psm[0] = P[0] # Reuse old value in case first data.. Psm[0] = P[0] # Reuse old value in case first data..
# ..point is skipped. Not plotted any way. # ..point is skipped. Not plotted any way.
P[-1] = P[-2] # Last data point (Nyquist f) can be invalid
# Loop through data points # Loop through data points
for x in range(x0, L): for x in range(x0, L):
@ -162,6 +162,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
fu = fNq # no data beyond fNq fu = fNq # no data beyond fNq
fl = fc**2 / fu # keep window symmetric fl = fc**2 / fu # keep window symmetric
# Find indices corresponding to frequencies
xl = bisect.bisect_left(freq, fl) # index corresponding to fl xl = bisect.bisect_left(freq, fl) # index corresponding to fl
xu = bisect.bisect_left(freq, fu) xu = bisect.bisect_left(freq, fu)
@ -178,7 +179,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
# Apply smoothing # Apply smoothing
Psm[x] = np.dot(g, P[xl:xu]) Psm[x] = np.dot(g, P[xl:xu])
if st == SmoothingType.levels: if st == SmoothingType.levels:
Psm = 10*np.log10(Psm) Psm = 10*np.log10(Psm)