0
$\begingroup$

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 enter image description here

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

PSD and RMS

$\endgroup$
0
$\begingroup$

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
    

`

$\endgroup$
7
  • $\begingroup$ Thanks. Just for clarity, here the T represents the number of time frames and not the total duration of the time. I cover 12-time frames covering a total duration of 2S. one more question, how can I obtain the single value RMS/power average in the frequency domain. Is it just the max of the power spectrum? $\endgroup$
    – Shan
    Apr 21 at 13:32
  • $\begingroup$ From the relation, after interchanging the last 2 lines, the Power_avg would be the PSD? $\endgroup$
    – Shan
    Apr 21 at 13:36
  • $\begingroup$ Yes, T represents the total number of time frames, I have edited to answer tor reflect that. As for obtaining a single value, the PSD is usually defined over the frequency axis. The RMS of the spectrum would be $\sqrt{\frac{1}{Nfft} \sum_{\omega=0}^{Nfft-1} P_x(\omega)^2}$, i.e, the square root of the mean of the squared values in $P_x$. $\endgroup$
    – orchi_d
    Apr 21 at 14:36
  • $\begingroup$ $Nfft$ is the total number of frequency bins, or size of the FFT $\endgroup$
    – orchi_d
    Apr 21 at 14:42
  • $\begingroup$ Thanks @Orchisama, I have added the RMS formula in the question and also changed the plot. PSD showing a Sine wave at 114dB 1Hz. I think the rms/power average should give me 114dB, but with the calculation i get 24dB. Can you please check, what i am missing here? $\endgroup$
    – Shan
    Apr 21 at 15:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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