2
$\begingroup$

So librosa.core.stft returns a complex single sided spectrogram.

My question is:

What normalization of the amplitude values should I perform afterwards?

I believe I have to multiply the amplitude outputs by 2 in order to preserve the energy that was assignated to the negative frequencies. I have not seen in the code that librosa.core.stft does this, although I might have missed it.

And then the other normalization step would be to divide the outputs by the length of the window size right?

However, after doing these two amplitude normalization steps, I still lack a factor of 2 somewhere. This can be seen in the picture below, where in the coloured picture the spectrogram has been calculated and an amplitude value of 0.5 is observed. And in the second plot, the FT of the whole signal has been computed and, applying the same normalization, the expected value of 1 is observed.

Picture: The frequency axis are not in Hz, the only important thing here are the amplitude values.

enter image description here

Here is the code:

%Import packages 
import scipy.fftpack as fft
import librosa 
import numpy as np

%Create sine wave of frequency 1000 Hz and 20000 samples 
frequency = 1000 
fs = 4 * frequency 
time_vector = 1/fs * np.arange(20000)

sine_wave = np.sin(2 * np.pi * frequency * time_vector)

window_size = 1024 
hop_length = 512

out  = librosa.core.spectrum.stft(sine_wave, n_fft = window_size, hop_length = hop_length)

%Take magnitude and normalize  
out = 2 * np.abs(out) / window_size

fig = plt.figure() 
axes = fig.add_subplot(111) 
im = axes.pcolormesh(out) 
fig.colorbar(im) 
plt.show()

%And plot also the fourier transform of the whole signal 
transform = 2 * np.abs(fft.fft(sine_wave)) / len(sine_wave)

%Plot only positive frequencies
plt.plot(transform[10000:])
$\endgroup$
  • $\begingroup$ I am not getting how to read the plot. An explanation would be helpful. Thankyou $\endgroup$ – Becky Sep 18 '18 at 19:36
  • $\begingroup$ @Becky : Welcome to SE.SP! Please do not add answers as comments. You will be able to comment once you have earned enough reputation. Please take some time to understand how SE.SP works, and enjoy! $\endgroup$ – Peter K. Sep 19 '18 at 11:35
5
$\begingroup$

I did answer a similar question a few years back, but can't find it. Basically, you are losing the energy because of the windowing. It's true to say that you should multiply the spectra by 2 in order to retrieve the energy from negative frequencies and divide it by the window size. However, this is only true for the rectangular window. The correct normalization should be following:

out = 2 * np.abs(out) / np.sum(window)

where np.sum(window) in case of rectangular (or no window in fact), is equal to the frame size. It kind of makes sense, right? If window attenuates the signal at its ends, then this energy should be recovered somehow.

Now back to your problem, librosa.stft by default uses hann window which looks like:

enter image description here

If you calculate the np.sum(np.hanning(window_size)) you will get $511.5$. If you read up more, the coherent gain of that window is approximately $0.5$. Which is defined as:

$$CG = \frac{\Sigma w(N)}{N} $$

So in order to solve your problem, I suggest to explicitly pass the window vector to the stft function:

import librosa
import numpy as np
import matplotlib.pyplot as plt

frequency = 1000
fs = 4 * frequency
time_vector = 1/fs * np.arange(20000)

sine_wave = np.sin(2 * np.pi * frequency * time_vector)

window_size = 1024
window = np.hanning(window_size)
hop_length = 512

out  = librosa.core.spectrum.stft(sine_wave, n_fft = window_size, hop_length = hop_length, window=window)
out = 2 * np.abs(out) / np.sum(window)

fig = plt.figure()
axes = fig.add_subplot(111)
im = axes.pcolormesh(out)
fig.colorbar(im)

plt.figure()
slice = out[:, 10]
plt.plot(slice)
plt.grid()

plt.show()

Which will return:

enter image description here

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.