a-weighting signal without digital filter

Question

I essentially want to a-weight a signal by using the fft in Python. But I guess I'm making some logical faults.

My Steps:

• Initialize Sine: y = np.array(np.rint(32768 * np.sin(1000 * 2.0 * np.pi * np.linspace(0,1,44100))))

• take RMS: print(20.0 *np.log10(rms_flat(y))) - gives me $87\textrm{ dB}$

• fft with hanning window:

  yf = np.abs(np.fft.fft(samples * np.hanning(len(samples))))

• take RMS: gives me around $130\textrm{ dB}$ - why is it not at $87\textrm{ dB}$ anymore?

• finally a-weight spectrum - see below - get Inf RMS...

import matplotlib.pyplot as plt
import numpy as np


#samplerate
Fs = 44100

#samplecount
A = Fs

#Sample Interval
Ts = 1/Fs

#Zeitvektor
t = np.linspace(0.0, Fs*Ts, Fs)

#Sinus Freq
ff = 1000

y = np.array(np.rint(32768 * np.sin(ff * 2.0 * np.pi * t)))

fig = plt.figure(figsize=(10,10))
fig.subplots_adjust(hspace=0.7)
fig.subplots_adjust(wspace=0.25)
plt.subplot(321)
plt.title("Samples")
plt.ylabel("Amplitude")
plt.xlabel("Time")
plt.plot(t, y)


def rms_flat(a):  # from matplotlib.mlab
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """
    return np.sqrt(np.mean(np.absolute(a)**2))


print(20.0 *np.log10(rms_flat(y)))


def fft(samples):
    yf = np.abs(np.fft.fft(samples * np.hanning(len(samples))))
    xf = np.linspace(0.0, 1 / (2.0 * Ts), Fs / 2)
    plt.subplot(322)
    plt.title("Abs Pos Spectrum")
    plt.ylabel("Amplitude")
    plt.xlabel("Frequency")
    plt.plot(xf, 2.0/A * yf[:A//2])
    print(20.0 * np.log10(rms_flat(yf)))
    return(yf, xf)

yf, xf = fft(y)


def aWeighting(yf, xf):
    f1 = 20
    f2 = 107
    f3 = 737
    f4 = 12194
    N = len(xf)
    i = 0
    yfa = []
    while i < N:
        a = 20.0 * np.log10((f4 ** 2 * xf[i] ** 4) / ((xf[i] ** 2 + f1 ** 2) * np.sqrt(xf[i] ** 2 + f2 ** 2) * np.sqrt(xf[i] ** 2 + f3 ** 2) * (xf[i] ** 2 + f4 ** 2))) + 2.0
        yfa.append(20.0 * np.log10(yf[i]) + a)
        i += 1

    yfa = np.array(yfa)
    print(rms_flat(yfa))
    plt.subplot(323)
    plt.title("A-weighted Spektrum in dB")
    plt.ylabel("dB")
    plt.xlabel("Frequency")
    plt.plot(xf, 2.0/A * yfa[:A//2])
    plt.show()
    return(yfa, weight)

yfa, weight = aWeighting(yf, xf)

Answer 1

The FFT has a processing gain of 10Log(N) that would apply to your measurements of the standard deviation, where N is the length of the signal, and the Hamming Window has a coherent gain of 0.54 (which is a loss). Refer to fred harris' famous paper On the Use of Windowing that details these gains and losses.

Let's see if this works out. You have 44100 samples, so the net processing gain is $10Log10(44100\times.54) = 43.8 \text{ dB}$

Your peak sinewave is 32768 so the rms should be $20log_{10}(32768/\sqrt{2})=87.3 \text { dB}$ which is consistent with your measurement.

After taking the window and the FFT, you saw the RMS signal increase to 130 dB, while you measured 87 dB before the FFT. 130-87 = 43 dB.

The difference of 43 dB is consistent with that predicted by the formula for the net processing gain given above.

---

The following is more details that may further help to understand how this occurs:

The coherent gain of your input signal is N. Your input signal in time is

$$A sin(\omega t) = \frac{A}{2}\frac{e^{j\omega t}}{j}-\frac{A}{2}\frac{e^{j\omega t}}{j}$$

The important item to see from Euler's expansion above is the sine is composed of two exponentials each with an amplitude of half of the sine wave itself. After the FFT, each of these is increased by a factor of N due to the coherent processing gain in the FFT.

Thus each of the two tones in the FFT would have a magnitude of $\frac{AN}{2}$ (if you sampled coherently at an exact integer number of samples this would be exact, otherwise we will see it slightly reduced from this due to "scalloping loss").

If sampled coherently, every other bin besides the two tones described will be zero (we could further show through Parseval's theorem the general case where all the energy in the signal is conserved, but to keep my explanation simple I will describe just the case of coherent sampling, meaning the sampling clock is an integer multiple of the sine wave).

Therefore when computing the rms, since there are only 2 non-zero values but we are using N samples, the rms calculation of the FFT would reduce to:

$$\sigma = \sqrt{\frac{\left(\frac{AN}{2}\right)^2+\left(\frac{AN}{2}\right)^2}{N}}$$

Which reduces to

$$\sigma = \sqrt{N}\frac{A}{\sqrt{2}}$$

Which we recognize as the $\sqrt{N}$ times the rms of our input sinusoid. When we window the signal before taking the FFT, the coherent amplitude (that previously increased by N) would now be reduced by the coherent gain of the window, which for the case of the Hamming window was 0.54.

Thus the standard deviation of the FFT is expected to be $\sqrt{(N\times0.54)}$ times larger than the standard deviation of the input sinusoid for the case of a Hamming Window being used. So in dB:

$$20log_{10}(\sigma) = 10log_{10}(N\times0.54) $$