user90664
user90664

Reputation: 75

Python histogram of RMS amplitude of audio file

I would like to make a histogram that bins RMS amplitude for an audio file. The goal is to show, for the full duration of the file, how much of the signal is at each amplitude level. Is there a Python package with a function for this? If not, how can it be coded?

I would also like to set the frequency range over which the analysis will be computed, e.g. between 1 and 6 kHz.

I have the following as a crude start, though I don't yet understand what it represents, and it certainly does not use RMS:

import numpy as np
import matplotlib.pyplot as plt
   
Fs, data = wavfile.read('file') 
print('data =',data)
print('number of samples in data =',len(data))

subset = data[0:44100] 
subset = abs(subset)
print('number of samples in subset =',len(subset))

plt.hist(subset, bins='auto')  
plt.show()

Upvotes: 0

Views: 1824

Answers (1)

Alexander Korovin
Alexander Korovin

Reputation: 1475

As far as I know, there is no special function in numpy for RMS, but you can do it like this

RMS = np.sqrt(np.mean(x**2))

And the question is, for which data (for which x) do you want to calculate the RMS. For example, you can apply RMS for each sample, than assuming you only have one channel in your wav-file:

length = data.shape[0] / Fs
print(f"length = {length}s")

RMS = lambda x: np.sqrt(np.mean(x**2))

sample = np.arange(int(length))
RMS_of_sample = np.zeros(sample.shape)
for ns in sample:
    # here you can apply the frequency window for the sample 
    RMS_of_sample[ns] = RMS(data[ns*Fs:(ns+1)*Fs])

plt.hist(RMS_of_sample, label="Left channel")
plt.show()

here you can also apply some signal windows. This code gives you something like this

enter image description here

for incoming signal:

enter image description here


ADDITION to the question in the comment regarding full/partial frequency range

If you want to analyze a complete signal in a certain frequency domain, you can apply, for example, simple filter (rectangular frequency window) for the frequency range [filter_freq_min, filter_freq_max] like this:

from scipy.fft import fft, ifft, fftfreq

filter_freq_min = 1000 # Hz
filter_freq_max = 2000 # Hz

freq = fftfreq(len(data), 1 / Fs)
data_fft = fft(data)

condition = np.logical_or(abs(freq) <= filter_freq_min, abs(freq) >= filter_freq_max)
filtered_data_fft = np.copy(data_fft)
filtered_data_fft [condition] = 0
filtered_data = np.real(ifft(filtered_data_fft ))

# show fft for incoming signal (blue) and filtered signal (orange) 
plt.plot(freq, np.abs(data_fft),'.')
plt.plot(freq, np.abs(filtered_data_fft ),'.')
plt.xlim( [10, Fs/2] )
plt.xlabel( 'Frequency (Hz)' )
plt.show()

# check RMS for filtered and unfiltered signal
print(RMS(filtered_data),RMS(data))

enter image description here

In this way, you can cycle through the required frequency ranges.

To play sound directly in Python, you can use

import sounddevice as sd # For playing/recording audio
sd.play(data, Fs)
sd.play(filtered_data, Fs)

Upvotes: 5

Related Questions