Блог

Я работаю над окончанием дипломной работы, в которой мне нужно измерить уровень звукового давления подводных записей (wav-файлов) на определенной частоте (2000 Гц). Итак, я придумал этот код:

«‘ def get_value(имя файла, 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)
 

»’
Однако я проверил значение с помощью audacity и оно совершенно другое.

Заранее благодарю.

Если вас интересует только одна частота, вам не нужно вычислять БПФ, вы можете просто использовать

 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)