1
$\begingroup$

I am trying to computing the overall spectrum of a given signal from its Short-Time-Fourier-Transform. I am using python to do this and I wrote my function to compute the STFT:

def compute_stft(x, window_size, hop_size):
# Calculate the number of frames == number of samples of the entire signal - framesize // hopsize + 1

num_frames = (len(x) - window_size) // hop_size + 1

# Initialize the STFT matrix [number of bins, number of frames]
stft_matrix = np.zeros((window_size // 2 + 1, num_frames))

# Generate the window function
window = sc.get_window("hann", window_size)

# Iterate over the signal in frames
for i in range(num_frames):
    start = i * hop_size
    end = start + window_size

    # Extract the current segment
    chunk = x[start:end]

    # Apply the window function
    windowed_segment = chunk * window
    # Compute the real-valued FFT of the segment
    stft_matrix[:, i] = np.abs(np.fft.rfft(windowed_segment)) * (2.0 / window.sum())
return stft_matrix

What I do next is to transpose the stft matrix for visualization reasons:

stft_matrix = compute_stft(x, 2**14, 512)
stft_matrix = stft_matrix.T
stft_freqs = np.linspace(0, fs / 2, stft_matrix.shape[1])

Then I would like to compute the overall spectrum by averagin between the time frames:

magnitude_spectrum = 20 * np.log10(stft_matrix.mean(axis=0) + 1e-60)

The I plot

plt.figure(figsize=(10, 6))
plt.semilogx(f_bins, magnitude_spectrum)
plt.xlim(100, fs / 2)
plt.ylim(-100, -6)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.title("Magnitude Spectrum")
plt.grid(True, which="both")

The result is the following enter image description here

While in Cool Edit Pro I obtain the following

enter image description here

And as you can see there is a very significant offset between the two graphs.. I suppose that it depends on how I am normalizing the STFT inside the function..but really I tried different kinds of normalization, even the one that is done inside scipy but without any considerable difference. Can someone help me understand How to solve this problem? Thank you in advance.

EDIT: $f_s = 48000 \texttt{Hz}$

$\endgroup$

2 Answers 2

4
$\begingroup$

You're almost there.

  • Your data should be scaled with dBFS reference (see the "reference" box in cool edit pro).

  • Cool Edit Pro probably computes the average power spectrum:

    1. Modify

      stft_matrix[:, i] = np.abs(np.fft.rfft(windowed_segment)) * (2.0 / window.sum()) 
      

      to

      stft_matrix[:, i] = 2.0 * ( np.abs(np.fft.rfft(windowed_segment)) / window.sum() ) ** 2 
      
    2. Take the mean and 10*np.log10():

      power_spectrum = 10 * np.log10(stft_matrix.mean(axis=0) + 1e-60)
      
$\endgroup$
5
  • $\begingroup$ Just to understand better, the normalization is just done by dividing with window.sum() and the you multiply by two everything in order to recorver the missing energy due to windowing? $\endgroup$ Aug 24, 2023 at 12:06
  • 1
    $\begingroup$ Yes, but don’t forget to square. The correct scaling is $$2 \cdot \left(\frac{|X|}{\sum w}\right)^2$$ $\endgroup$
    – Jdip
    Aug 24, 2023 at 12:08
  • $\begingroup$ And you had $$2 \cdot \frac{|X|}{\sum w}$$ $\endgroup$
    – Jdip
    Aug 24, 2023 at 12:11
  • $\begingroup$ Of course! As far as you know, is this normalization valid for any kind of window? Also, If I wanted to plot the magnitude spectrum (i.e. not to be compared with Cool Edit Pro) would it be that way I proceeded originally correct? $\endgroup$ Aug 24, 2023 at 12:16
  • $\begingroup$ There are subtleties outside the scope of this question, but in general, yes to both these questions ;) $\endgroup$
    – Jdip
    Aug 24, 2023 at 12:21
1
$\begingroup$

A few things here

  1. The overall spectrum is typically determine by averaging in energy, not in amplitude.
  2. When you plot something in dB you need to be clear about what reference you are using. That depends a bit on the units of your original data. Could be dBV, dBSPL, dBFS, normalized etc. This will determine the offset. It also depends on the DFT scaling convention used. All of this will determine the absolute offset of the graph. Cool edit is most likely to us dBFS (0dB = digital full scale).
$\endgroup$

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.