STFT computation

I'm trying to compute the STFT by my self. Thus, I wrote the following code:

hann_win = scipy.signal.hanning(256, sym=True)
frame_spectrum = np.zeros((int(fft_Size / 2 + 1), N_frames))
for i in range(N_frames):
frame_spectrum[:, i] = abs(np.fft.rfft(hann_win * signal[start:start + frame_size], fft_Size))**2
start += 128


fft_size = frame_size = 256. Using the scipy.signal.stft function, the output result is different. It looks like it is normalized in some way.

f, t, STFT = scipy.signal.stft(signal, Fs, 'hann', 256, 128, 256, boundary = None, padded = False, axis = 0)


The ouput for my implementation is:

frame_spectrum
Out[239]:
array([[1.39350590e+05, 1.28392811e+05, 3.01856941e+05, ...,
6.44869243e+05, 6.86848227e+05, 2.35430292e+05],
[2.45825457e+04, 1.11256638e+04, 5.85035486e+04, ...,
2.01218021e+05, 1.96921261e+05, 1.20836787e+05],
[2.12379247e+04, 2.13193995e+04, 6.98189265e+03, ...,
4.81385299e+03, 2.79953520e+04, 1.11993281e+04],
...,
[2.35902280e+02, 2.06442959e+02, 8.19946703e+02, ...,
5.85829771e+02, 4.35596749e+01, 3.36198896e+02],
[2.77258143e+02, 1.61436920e+02, 6.49894194e+02, ...,
5.80593478e+02, 3.19306802e+02, 1.09445294e+03],
[3.93939998e+02, 1.50623962e+02, 5.51785739e+02, ...,
4.90244067e+02, 4.83394024e+02, 1.97937846e+03]])


For the scipy.signal.stft function, the output is:

STFT
Out[242]:
array([[ 2.9267783 +0.        j,  2.8191867 +0.        j,
4.3128223 +0.        j, ...,  6.2927547 +0.        j,
6.5078945 +0.        j,  3.778513  +0.        j],
[-1.0943314 -0.51656187j, -0.8138022 +0.10491034j,
-1.8788935 -0.0876543 j, ..., -3.498077  -0.11071016j,
-3.3105822 +1.0691631 j, -1.8662903 -1.960599  j],
[-0.8903561 +0.71827334j, -1.1007539 +0.3362887 j,
-0.60977626+0.25182465j, ...,  0.5338975 +0.02365825j,
-0.07851043-1.3080015 j,  0.15810308+0.80889297j],
...,
[-0.06972471+0.09730107j, -0.06172764-0.09533134j,
-0.22229262+0.02101982j, ..., -0.18767883+0.02727323j,
-0.05091092-0.00271361j, -0.12367946-0.07388743j],
[ 0.12895514-0.01280107j,  0.09401653-0.03171136j,
0.19183227+0.05289177j, ...,  0.18554814-0.02986654j,
0.12113474+0.06995475j,  0.2588121 -0.00447672j],
[-0.15553662+0.        j, -0.09567829+0.        j,
-0.18430477+0.        j, ..., -0.17289034+0.        j,
-0.17330773+0.        j, -0.34770313+0.        j]],
dtype=complex64)


Could anyone give me an explanation please? I missed something?

Edit:

abs(STFT)**2
Out[3]:
array([[8.56603146e+00, 7.94781351e+00, 1.86004372e+01, ...,
3.95987625e+01, 4.23526917e+01, 1.42771597e+01],
[1.46439719e+00, 6.73280120e-01, 3.53792381e+00, ...,
1.22487984e+01, 1.21030636e+01, 7.32698822e+00],
[1.30865073e+00, 1.32474923e+00, 4.35242742e-01, ...,
2.85606265e-01, 1.71703196e+00, 6.79304421e-01],
...,
[1.43290339e-02, 1.28983660e-02, 4.98558395e-02, ...,
3.59671712e-02, 2.59928545e-03, 2.07559634e-02],
[1.67932957e-02, 9.84471757e-03, 3.95971648e-02, ...,
3.53201218e-02, 1.95672903e-02, 6.70037419e-02],
[2.41916403e-02, 9.15433560e-03, 3.39682512e-02, ...,
2.98910681e-02, 3.00355703e-02, 1.20897464e-01]], dtype=float32)


Ratio between frame_spectrum and abs(STFT)**2:

absolute_frame_spectrum / (abs(STFT)**2)
Out[59]:
array([[16267.81209577, 16154.4821358 , 16228.4863926 , ...,
16285.0857511 , 16217.34535357, 16489.99500715],
[16786.8020194 , 16524.568988  , 16536.12449353, ...,
16427.57227861, 16270.36492832, 16492.01333113],
[16228.87161575, 16093.15862715, 16041.37639628, ...,
16854.8578121 , 16304.50260254, 16486.46436409],
...,
[16463.23688692, 16005.35752565, 16446.35236183, ...,
16287.90231531, 16758.32674005, 16197.70134645],
[16510.04947047, 16398.32925761, 16412.64461434, ...,
16438.03727133, 16318.39649219, 16334.20625631],
[16284.13753068, 16453.83875639, 16244.16092896, ...,
16401.02205314, 16094.05180593, 16372.37367461]])


Spectrogram for my implementation: Spectrogram for scipy.signal.stft function:

So, after a search in the official sources of the scipy.signal.stft, I found the normalization.

They compute the scale factor as following:

scale = 1.0 / win.sum()**2


where 'win' represents the window function selected. Then, if the 'stft' mode is selected, the root from the scale is extracted:

 scale = np.sqrt(scale)


Now, I tried in the same way, and the results are almost the same. My implementation:

hann_win = scipy.signal.hamming(256, sym=True)
for i in range(N_frames):
frame_Spectrum[:, i] = np.fft.rfft(hann_win * signal[start:start + frame_size], fft_Size)
scale = 1.0 / ((hann_win).sum() ** 2)
scale = np.sqrt(scale)
frame_Spectrum *= scale
frame_Spectrogram[:, i] = abs(frame_Spectrum[:, i])**2
step += 128


The scipy.signal.stft version:

f, t, STFT = scipy.signal.stft(signal, Fs, 'hamm', 256, 128, 256, boundary = None, padded = False, axis = 0)


The output for both versions is listed below. For my implementation:

frame_Spectrogram
Out[28]:
array([[1.11538752e+03, 8.07521488e+01, 1.48759566e+02, ...,
1.42030962e+02, 8.61416251e+01, 1.29774958e+02],
[4.34235618e+02, 4.57448938e+01, 6.62891618e+01, ...,
1.34617117e+01, 7.47917194e+01, 1.41128514e+02],
[6.54247473e+01, 2.09806259e+01, 5.72439770e+00, ...,
2.72848404e+01, 9.29473602e+01, 9.99377515e+01],
...,
[2.34281203e-02, 2.17043722e-02, 5.95170905e-01, ...,
4.17076387e-02, 2.09002175e-01, 1.30255452e+00],
[5.13139103e-02, 9.10540408e-02, 3.43592577e-02, ...,
2.47531581e-02, 7.68725761e-02, 3.16595192e-01],
[5.54714753e-04, 1.41252517e-01, 3.66720655e-02, ...,
1.23587668e-01, 1.65866874e-02, 3.71481824e-02]])


For the scipy.signal.stft function:

abs(STFT)**2
Out[29]:
array([[1.11314111e+03, 7.97755508e+01, 1.49654175e+02, ...,
1.41768829e+02, 8.52271576e+01, 1.29163040e+02],
[4.31478424e+02, 4.48896141e+01, 6.66922226e+01, ...,
1.31336327e+01, 7.38092880e+01, 1.40289795e+02],
[6.42291183e+01, 2.08642082e+01, 5.85985231e+00, ...,
2.68950195e+01, 9.29605713e+01, 9.92776413e+01],
...,
[2.31771655e-02, 1.91425439e-02, 5.86372852e-01, ...,
3.92387286e-02, 2.08843708e-01, 1.30342722e+00],
[5.06063662e-02, 8.90289247e-02, 3.40913758e-02, ...,
2.56771799e-02, 7.78447017e-02, 3.19540411e-01],
[6.48763496e-04, 1.38118595e-01, 3.60494740e-02, ...,
1.23369396e-01, 1.58909075e-02, 3.86624187e-02]], dtype=float32)


As you can see, the results have the same order now and and they differ by several decimals.

The $\texttt{scipy.signal.stft()}$ function returns the complex-valued STFT. You can see in the output that its datatype is $\texttt{complex64}$. You however, are taking the square of the absolute value of the complex-valued STFT coefficients, which is typically called the spectrogram.

• 𝚜𝚌𝚒𝚙𝚢.𝚜𝚒𝚐𝚗𝚊𝚕.𝚜𝚝𝚏𝚝() function returns the complex values of STFT, it is correct that the results shown by me in the second case are not the square of the absolute values. Even so, if you compute the square of the absolute values listed above, you will see that are not the same with the values resulted from my implementation. Feb 21 '18 at 6:05
• It would be helpful to visualise the results, ie. by plotting the two spectrograms side by side. It seems your own values are quite a bit larger than the values given by the function. You could try to compute the ratio between your own values and the other ones, and see if it constant. Then it is just a matter of scaling. It might be that the scipy function scales the analysis window such that its total energy or the sum of its absolute values is 1. Or it might be scaled such that when you sum all the time-shifted analysis windows, they sum to 1 at all sample positions. Feb 21 '18 at 9:27
• I already did the the ratio between the square of the absolute values of the frame spectrum and the square of the absolute values resulted from the STFT function. I have seen that the ratio is almost constant, but I didn't found an explanation for how the normalization factor has been choose. Also I analyzed the spectrogram, and both look pretty much the same. On the scipy.signal.stft official documentation they didn't say about any normalization. Check my last edit to view the spectrograms and the discussed ratio. Thank you! Feb 21 '18 at 9:44