如何获得声音文件特定频率的功率？

Question

我的学位论文即将结束，我必须在特定频率 (2000Hz) 下测量水下录音（wav 文件）的声压级。所以我想出了这个代码：

''' def get_value(filename, f0, NFFT=8192, plot = False):

#Load audio
data, sampling_frequency = soundfile.read(filename)

# remove stereo
if len(data.shape)> 1:
    data = data[:, 0]

# remove extra length
if len(data)>sampling_frequency:
    data = data[0:sampling_frequency]

# remove DC
data = data - data.mean()

# power without filtering
total_power = 10*np.log10(np.mean(data**2))
    
# fft
NFFT = 4096 # number of samples in the FFT
window = np.array(1) #np.hamming(len(data))
fftdata = np.fft.fft(data / NFFT, n = NFFT)

SPL = 20 * np.log10(np.abs(fftdata))            # Sound Pressure Level [dB]
freq = np.linspace(0, sampling_frequency, NFFT) # frequency axis [Hz]

# take value at desired frequency
power_at_frequency = SPL[np.argmin(np.abs(freq-f0))]
print(power_at_frequency)

''' 但是，我大胆地检查了这个值，结果完全不同。

先谢谢了。

Answer 1

如果您只对一个频率感兴趣，则无需计算 FFT，您可以简单地使用

totalEnergy = np.sum((data - np.mean(data)) ** 2)
freqEnergy  = np.abs(np.sum(data * np.exp(2j * np.pi * np.arange(len(data)) * target_freq / sampling_freq)))

如果您使用 FFT 并且窗口大小不是波周期的倍数，则频率将泄漏到其他频率。为避免这种情况，您的

import numpy as np;
import matplotlib.pyplot as plt
sampling_frequency = 48000;
target_frequency = 2000.0;

ns = 1000000;

data = np.sin(2*np.pi * np.arange(ns) * target_frequency / sampling_frequency);

# power
print('a sine wave have power 0.5 ~', np.mean(data**2), 'that will be split in two ')

## Properly scaled frequency
plt.figure(figsize=(12, 5))
plt.subplot(121);
z = np.abs(np.fft.fft(data[:8192])**2) / 8192**2
print('tuned with 8192 samples', max(z), ' some power leaked in other frequencies')
plt.semilogy(np.fft.fftfreq(len(z)) * sampling_frequency, z)
plt.ylabel('power')
plt.title('some power leaked')

plt.subplot(122);
# 6000 samples = 1/8 second is multiple of 1/2000 second
z = np.abs(np.fft.fft(data[:6000])**2) / 6000**2
print('tuned with 6000 samples', max(z))
plt.semilogy(np.fft.fftfreq(len(z)) * sampling_frequency, z)
plt.xlabel('frequency')
plt.title('all power in exact two symmetric bins')
## FFT of size not multiple of 2000
print(np.sum(np.abs(np.fft.fft(data[:8192]))**2) / 8192)