Skip to content

Commit

Permalink
resolve spectrum estimation 'wrong' #947
Browse files Browse the repository at this point in the history 
scale DC and Nyquist frequency by just /N not 2/N
  • Loading branch information
@sammlapp
sammlapp committed May 2, 2024
1 parent 8f383df commit a9cdf88
Showing 1 changed file with 10 additions and 7 deletions.
17 changes: 10 additions & 7 deletions opensoundscape/audio.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -732,17 +732,20 @@ def spectrum(self):

# Compute the fft (fast fourier transform) of the selected clip
N = len(self.samples)
fft = scipy.fft.fft(self.samples)
fft = scipy.fft.rfft(self.samples)
fft = np.abs(fft) # get the magnitude of the fft

# create the frequencies corresponding to fft bins
freq = scipy.fft.fftfreq(N, d=1 / self.sample_rate)
freq = scipy.fft.rfftfreq(N, d=1 / self.sample_rate)

# remove negative frequencies and scale magnitude by 2.0/N:
fft = 2.0 / N * fft[0 : int(N / 2)]
frequencies = freq[0 : int(N / 2)]
fft = np.abs(fft)
# scale magnitude by 2.0/N,
# except for the DC and sr/2 (Nyquist frequency) components
fft *= 2.0 / N
fft[0] *= 0.5
if N % 2 == 0:
fft[-1] *= 0.5

return fft, frequencies
return fft, freq

def normalize(self, peak_level=None, peak_dBFS=None):
"""Return audio object with normalized waveform
Expand Down

1 comment on commit a9cdf88

@WMXZ-EU
Copy link

@WMXZ-EU WMXZ-EU commented on a9cdf88 May 3, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I understand now, and I may have formalated the issue wrong: It is the energy that is maintained (spectral power) and not the amplitude
I did the folloing test

import numpy as np
import scipy.fft
x=np.array([0,0,4,0,0,0,0])
N=len(x)

def Power(x,N): return np.sum(np.abs(x)**2,axis=0)/N

x1=scipy.fft.fft(x)
print(N,x1)
print(Power(x,1),Power(x1,N))

x2=scipy.fft.rfft(x)

n2=N//2
if N%2==1: n2 +=1

x2[1:n2] *=np.sqrt(2)
print(N//2,x2)
print(Power(x,1),Power(x2,N))

to demonstrate the proper scaling

In your case I would suggest
N = len(self.samples)
fft = scipy.fft.rfft(self.samples)

fft = np.abs(fft) **2/ N # scale here all spectral bins

n2=N//2
if N%2==1: n2 +=1 # estimate number of frequencies excluding Nyquist (happens when N is even)

fft[1:n2] *=2 # correct for negative frequencies

and Level_dB = 10*np.log10(fft)

Please sign in to comment.