How can I create a histogram showing the RMS amplitude for each frequency bin of an audio file

Using Python or Matlab, how can I create a histogram in which each bin is equal to a proportional frequency range (e.g. one octave) and the data added to each bin is a count of the number of samples that have any amplitude within that frequency range?

Here's what I have so far:

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:
RMS_of_sample[ns] = RMS(data[ns*Fs:(ns+1)*Fs])

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


But this seems to create a bin equal to the full spectral range, thus only one bar appears in the histogram. I'm also not sure that the data in the bar actually represent counts of amplitude occurrences.

If what you want is to determine the amount of energy your signal has in particular frequency bands you can perform a FFT of your signal. For that you need to make sure that you have amplitude data and not energy data. So you would want to have sound pressure Pa and not Pa^2. After having calculated a FFT with that the values have to be multiplied by $$1/\sqrt{2}$$ to get to the RMS-format.

Depending on your time data length and sampling time you get a specific frequency resolution of your signal in the frequency domain. With $$N$$ time samples and a sampling time of $$\Delta t$$ you get a frequency resolution of $$1/(\Delta t \cdot N)$$ Herz. Now if you want to know the energy contribution of a specific frequency range you can calculate the RMS-value of that frequency range using the earlier calculated linear spectra in RMS-format. According to Parseval's theorem one can use the following equation to calculate the RMS-value of a specific frequency range: Therein the $$A_k$$s are the amplitudes of different frequencies in the spectra:

I used that approach to calculate the Energy of 10 Hz bins of a signal containing two sines with amplitude $$2$$ and Frequencies $$5$$ and $$45$$ Hz. Sampling frequency is 0.0098.

Here is the time Signal:

The resulting linear Spectra with Peak to Peak values is:

From that I calculated the RMS-Values of $$\Delta 10$$ Hz bins. So the first bar indicates the RMS-Value of the signal from $$0 - 10$$ Hz. The second from $$10 - 20$$ Hz and so on.

If some noise is added the spectra changes to:

The RMS-values of the $$\Delta 10$$ Hz bins become then:

I further suggest looking at Root Mean Square and Overall Level

Here is the Matlab Code I used to do the above (probably still containing errors though):

clear all
close all

t=linspace(0,10,1024); % time vector
y=2*sin(5*2*pi*t)+2*sin(45*2*pi*t)+randn(1,length(t)); % time signal
Fs=1/(t(2)-t(1));
Fn=Fs/2;
N=length(y);
yf=fft(y,N);
df = Fs/N; % Frequency Resolution

plot(t,y)
grid on
title('Time Signal')
ylabel('Amplitude')
xlabel(['Time Resolution: ',num2str(1/Fs),' s'])

amplH = abs(yf);
amplitudengang = fftshift(amplH/N);

x_fn = 0 : df : Fn-df;
x_fa = 0 : df : Fs-df;

figure
stem(x_fa-Fn, amplitudengang, 'b.-')
axis([-Fn Fn 0 max(amplitudengang)])
title('Two sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),' Hz'])
grid

amplitudengang=[amplitudengang(513) amplitudengang(514:end).*2];

figure
stem([0:df:(Fn-df)], amplitudengang, 'b.-')
axis([0 Fn 0 max(amplitudengang)])
title('Single Sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),'Hz'])
grid

amp10=amplitudengang(1:round((10-df)/df));
amp20=amplitudengang(round((10-df)/df):round((20-df)/df));
amp30=amplitudengang(round((20-df)/df):round((30-df)/df));
amp40=amplitudengang(round((30-df)/df):round((40-df)/df));
amp50=amplitudengang(round((40-df)/df):round((50-df)/df));

rms10=amp10./sqrt(2);
rms20=amp20./sqrt(2);
rms30=amp30./sqrt(2);
rms40=amp40./sqrt(2);
rms50=amp50./sqrt(2);

rms10=sqrt((rms10(1)/2)^2+sum(rms10(2:(end-1)).^2)+(rms10(end)/2)^2);
rms20=sqrt((rms20(1)/2)^2+sum(rms20(2:(end-1)).^2)+(rms20(end)/2)^2);
rms30=sqrt((rms30(1)/2)^2+sum(rms30(2:(end-1)).^2)+(rms30(end)/2)^2);
rms40=sqrt((rms40(1)/2)^2+sum(rms40(2:(end-1)).^2)+(rms40(end)/2)^2);
rms50=sqrt((rms50(1)/2)^2+sum(rms50(2:(end-1)).^2)+(rms50(end)/2)^2);

figure
stem([5 15 25 35 45],[rms10,rms20,rms30,rms40,rms50])
axis([0 Fn 0 max([rms10,rms20,rms30,rms40,rms50])])
title('RMS values per 10Hz')
ylabel('Amplitude')
xlabel(['Frequency in Hz'])
grid on


Edit

To reduce leakage the time data needs to be windowed. After that the resulting spectra needs to be multiplied by the matching correction factor:

• Without proper windowing in the time domain, this method can create significant spectral leakage, i.e. a sine wave at 1 kHz can create non trivial energy in, say, the 250 Hz octave bin. – Hilmar Mar 8 at 13:10
• That's true indeed. In my example in the spectra the sine at 45 Hz is significantly more spread than the one at 5 Hz. As you mentioned that can be accounted for with windowing. So the time data has to be treated with a window and then the spectra has to be multiplied with a correction factor. – Matthias La Mar 8 at 13:16

There are two different ways to do this:

1. Bandpass filter and then calculate RMS
2. Do short term Fourier transform and

Unless you have a really large number of bands, method number 1 is the better one. Find yourself some octave band filters online or design them yourself. There is one filter for each frequency band. To get the energy in the 1 kHz band simply filter the input with a 1 kHz octave band filter.

Then apply a "running" RMS to the output. There are lot of different ways to tune the time domain behavior. The simplest would be to chop it in frames and calculate the RMS in each frame. You can also do more fancy stuff: a very quick attack and a slow decay, which is would a typical "meter bridge" would do.

For the FFT method you would chop the signal into overlapping frames, apply a window, do an FFT. Then apply a different spectral mask for each octave (or 3rd octave, etc) and calculate the (complex) RMS in each frequency band. The spectral masks should be derived by looking at ANSI S1.11-2004 , which is a standard for "acceptable" filters. This method has less flexibility in adjusting the time domain properties of the output, you simply get one RMS value per band and per frame.