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)