Overall spectrum computation from Short-Time-Fourier-Transform

Question

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

While in Cool Edit Pro I obtain the following

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}$

Answer 1

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)
     

Answer 2

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).