Calculate average from a power spectrum

Question

I want to calculate the average of a power spectrum in the complete frequency range. I have done a stft and want to calculate power in dB and then the average. Please let me know from the code below if my understanding is correct and how to calculate the average of the power spectrum.

y, sr = librosa.load('Noise Recorder.wav',sr=None,duration=2) 
n_fft = 8192
S = librosa.stft(y=y,n_fft=n_fft,win_length=n_fft,hop_length=n_fft,window='hann')     # Short FFT
F = librosa.fft_frequencies(sr=sr, n_fft=n_fft) # Frequency components
Abs = np.abs(S)  # Same as mag, phase = librosa.magphase(S)
Abs_avg = np.mean(Abs, axis=1)
Power_avg =librosa.power_to_db(Abs_mean**2) # dB
Scale_dB = 86.85    # Scaling factor

Power_Scaled_dB = Power_avg + Scale_dB  #Plot #Max:114dB at 1KHz
RMS = np.sqrt(np.sum(Power_Scaled_dB**2) / n_fft) # n_fft=8192 #Ans:24.4

###### Additional RMS tryouts #####
# Method 1 : Based on PSD formula
Px = 10**(np.max(114/10)) # Ans: 251188643150.9
RMS_Px = np.sqrt(np.sum(Px**2) / n_fft) # Ans: 2775268639.5
RMS_Px_dB = 10*np.log10(RMS_Px) # Ans: 94

# Method 2: Based on amplitude fomrula
A = np.sqrt(Px) # Ans: 501187.2
RMS_A = A / np.sqrt(n_fft) # nfft= 8192(fft length), Ans: 5537.388930343294
RMS_A_dB = 20*np.log10(RMS_A) # Ans: 74

With a 114dB 1Hz Sine wave input signal, after stft i obtain the psd in the plot

Now I want to derive the PSD from the spectrum, to further calculate the average/RMS of the spectrum. Given by the relation PSD and from there I want to calculate RMS given by the relation RMS/Leq. I am unsure how to derive the PSD from the Power_avg

and

Answer 1

By Welch's method, you can calculate the power spectrum by averaging the magnitude of a bunch of FFT frames (or time frames of the STFT). If your signal is $x(t)$, and its STFT is $X(\omega, \tau)$, where $\omega$ is the frequency bin, $\tau$ is the time frame and $T$ is the total number of frames, the PSD is

$$P_x(\omega) = \frac{1}{T} \sum_{\tau=0}^{T-1} |X(\omega, \tau)|^2 $$

When calculating the STFT, you need to make sure you get constant overlap add. This means, with the Hann window, you need a hopSize of N/2 (where N is your window length). I'd also keep the window length shorter, about 256/512 samples. This is not compulsory in your case, because your signal is stationary. Had it not been stationary you should have chosen a smaller window length. Here is your code, modified.

    y, sr = librosa.load('Noise Recorder.wav',sr=None,duration=2) 
    n_fft = 8192
    win_len = 512
    S = librosa.stft(y=y,n_fft=n_fft,win_length=win_len,hop_length=win_len/2,window='hann')     # Short FFT
    F = librosa.fft_frequencies(sr=sr, n_fft=n_fft) # Frequency components
    Abs = np.abs(S)  # Same as mag, phase = librosa.magphase(S)
    
    PSD = np.mean(Abs**2, axis=1).  #this is the power spectrum
    PSD_dB = 10*np.log10(PSD). #power spectrum in dB
    RMS = np.sqrt(np.mean(PSD**2))) #RMS of the power spectrum
    RMS_dB = 20*np.log10(RMS) # RMS of power spectrum in dB

    ##Theoretical RMS (in dB) for a sine wave of amplitude A dB
    A_lin = 10**(A/20.0)           #convert dB to linear amplitude, A = 114 in your case
    RMS_exp = 20*np.log10(A_lin/n_fft)   #expected RMS in dB
    

`