0
$\begingroup$

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 my implementation: Spectrogram for scipy.signal.stft function: Spectrogram for scipy.signal.stft function:

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 𝚜𝚌𝚒𝚙𝚢.𝚜𝚒𝚐𝚗𝚊𝚕.𝚜𝚝𝚏𝚝() 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. $\endgroup$
    – Alexandru Rusu
    Commented Feb 21, 2018 at 6:05
  • $\begingroup$ 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. $\endgroup$
    – hulappa
    Commented Feb 21, 2018 at 9:27
  • $\begingroup$ 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! $\endgroup$
    – Alexandru Rusu
    Commented Feb 21, 2018 at 9:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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