|
|
"""
This Code can calculate the Discrete Fourier Transformation (DFT) of a Noise Signal of different formats calculate theRMS Noise level over a given bandwidth bandwidth.
__author__ = "Lenny Hoenen"__email__ = "l.hoenen@posteo.de", "lennart.hoenen@stud.uni-heidelberg.de"
"""
import numpy as npimport matplotlib.pyplot as pltfrom scipy.fft import rfft, rfftfreqfrom scipy.integrate import simpsimport scipy.io.wavfile as wavfileimport json
def Average_Fourier_Uniform(Data, Samplingrate): """
Calculate the average FFT of multiple real, one dimensional signals with a given samplingrate and samplesize N. Returns the frequency range xf and the Linear Spectral Density of the signal over that range yf. Data is of the shape M x N, where M is the number of measurements to average and N is the Number of Data points in each measurement. The designation "Uniform" indicates, that no window function is used for the DFT. For Time Series consisting purely of noise, this should give the optimal results. """
N = len(Data[0]) T = 1/Samplingrate yf_av = [] for i in Data: yf = rfft(i) yf = (2 * np.abs(yf)**2) / (Samplingrate * N) #Calculates the Power Spectral Density of one measurement/signal see Lenny's Thesis 3.2.4 yf_av.append(yf) yf_av = np.sqrt(np.average(yf_av, axis=0)) #Average and convert from Power Spectral Density to Linear Spectral Density see Lenny's Thesis 3.2.4 xf = rfftfreq(N, T) return xf, yf_av
def Calc_Power_Noise_RMS(Frequencies, LSD, maxFreq): """
Calculate the RMS Noise Power over a specified frequency range. As the impedance of the coils acts as a low pass filter, the current and thus the magnetic field cannot follow arbitrarily large frequencies. The Osci might measure high freq noise in the monitoring voltage, but this will hardly exist in the magnetic field. Therefore, one might want to limit the bandwidth to calculate the RMS noise in the current and thus in the magnetic field. The below functions allow to select a maximum Frequency (maxFreq) and calculate the RMS noise accordingly. LSD is the Linear Spectral Density as given by Average_Fourier_Uniform. For Discussions of the stability of a magnetic field, the Voltage Noise Density and its RMS value will be most meaningful. The result of the function Calc_Voltage_noise_RMS() is therefore more meaningful. See Lenny's Thesis 3.4
The desired max. frequency might not be part of the list of discrete frequencies in the space of the DFT. Therefore a closest match is found and selected. """
res = max(Frequencies) index = len(Frequencies)
for i in range(len(Frequencies)): if res > np.abs(maxFreq - Frequencies[i]): res = np.abs(maxFreq - Frequencies[i]) index = i
print("You chose the maximum Frequency ", maxFreq, "Hz. The closest matching frequency contained in the DFT is ", Frequencies[index], "Hz with index ", index, ".") PowerNoiseRMS = simps(LSD[0:index]**2, Frequencies[0:index])
return PowerNoiseRMS
def Calc_Voltage_Noise_RMS(Frequencies, LSD, maxFreq): """
See explanation of Calc_Power_Noise_RMS() and Lenny's Thesis 3.4. """
PowerNoiseRMS = Calc_Power_Noise_RMS(Frequencies, LSD, maxFreq) VoltageNoiseRMS = np.sqrt(PowerNoiseRMS)
return VoltageNoiseRMS
"""
HERE IS JUST SOME DIFFERENT WAYS TO OPEN OSCI DATA:
#OPEN wav files (wav is quite efficient and easy to handle for noise analysis)Samplingrate, Data = wavfile.read("Z:/Directory/File.wav")Data = Data * AbsoluteVoltageRange / NumberofBitsWith the Handyscope HS6Diff, the minimum Voltage Range was +/-200mV. Therefore the AbsoluteVoltageRange is 400mV. The highest resolution is 16Bit. Therefore the NumberofBits is 2**16
#OPEN JSON filesf = open("Z:/Directory/File.json")jdata = json.load(f)[0]joutput = jdata["outputs"]JDATA = joutput[0]["data"]f.close()JDATA = np.array(JDATA)
#OPEN CSV filesdata = np.genfromtxt("Z:/Directory/file.csv", skip_header=9)data = data.T"""
"""
EXAMPLE:"""
samplingrateFinal, dataFinal = wavfile.read("Z:/Users/Lenny/Power Supply Mk.2/Noise Measurements/NoiseDensityHHCoil2ndPatchedChannelBatteryInput.wav")dataFinal = dataFinal*0.4/2**16 #Correct conversion from bit values to physical data, for the used Voltae range and precision.dataFinalshaped = np.reshape(dataFinal, [5,4000000]) #Shape one long continuous measurement into 5 shorter measurements to be averaged
samplingrateprototype, dataprototype = wavfile.read("Z:/Users/Lenny/Power Supply Mk.2/Noise Measurements/NoiseDensityHHCoilPrototypeBatteryInput.wav")dataprototype = dataprototype*0.4/2**16dataprototypeshaped = np.reshape(dataprototype, [5,4000000])
samplingratebackground, databackground = wavfile.read("Z:/Users/Lenny/Power Supply Mk.2/Noise Measurements/NoiseDensityHHCoilsFinalPowerSupplyWithPIBatteryInputBUTEVERYTHINGTURNEDOFF.wav")databackground = databackground*0.4/2**16databackgroundshaped = np.reshape(databackground, [5,4000000])
x_Final_av, y_Final_av = Average_Fourier_Uniform(dataFinalshaped, samplingrateFinal)x_prototype_av, y_prototype_av = Average_Fourier_Uniform(dataprototypeshaped, samplingrateprototype)x_background_av, y_background_av = Average_Fourier_Uniform(databackgroundshaped, samplingratebackground)print(Calc_Voltage_Noise_RMS(x_Final_av, y_Final_av, 30000))print(Calc_Voltage_Noise_RMS(x_prototype_av, y_prototype_av, 30000))print(Calc_Voltage_Noise_RMS(x_background_av, y_background_av, 30000))
plt.rcParams["figure.figsize"] = [7,5]plt.rcParams["font.size"] = 12plt.figure(2, figsize=[12,7])plt.loglog(x_background_av, y_background_av, linewidth=1, label="Background Noise in Lab", alpha=0.7)plt.loglog(x_prototype_av, y_prototype_av, linewidth=1, label="Noise Prototype in carton box", alpha=0.7)plt.loglog(x_Final_av, y_Final_av, linewidth=1, label="Noise Final Powersupply in carton box", alpha=0.7)#plt.loglog(x_DC1V_av, y_DC1V_av, linewidth=1, label="Noise with FreqGen and Shielding", alpha=0.7)plt.legend()plt.grid(which="major")plt.grid(which="minor", alpha=0.2)plt.xlabel("Frequency [Hz]")plt.ylabel(r"Voltage Noise Spectrum [$V/ \sqrt{Hz}$]")plt.title("5 Measurements of 4M Samples at 6.52MHz averaged")plt.show()plt.close()
|
