Commented and Cleaned version of the Noise Analysis Code that was used in Lenny's Thesis.

FourierTrafoAndRMSNoiseAnalysis.py

@@ -0,0 +1,136 @@
+"""
+This Code can calculate the Discrete Fourier Transformation (DFT) of a Noise Signal of different formats calculate the
+RMS Noise level over a given bandwidth bandwidth.
+
+
+__author__ = "Lenny Hoenen"
+__email__ = "l.hoenen@posteo.de", "lennart.hoenen@stud.uni-heidelberg.de"
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fft import rfft, rfftfreq
+from scipy.signal.windows import blackman
+from scipy.integrate import simps
+import scipy.io.wavfile as wavfile
+import 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(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(LSD)
+    index = len(LSD)
+
+    for i in range(len(LSD)):
+        if res > np.abs(maxFreq - LSD[i]):
+            res = np.abs(maxFreq - LSD[i])
+            index = i
+
+    print("You chose the maximum Frequency ", maxFreq, "Hz. The closest matching frequency contained in the DFT is ", LSD[index], "Hz with index ", index, ".")
+    PowerNoiseRMS = simps(y_Final_av[0:index]**2, x_Final_av[0:index])
+
+    return PowerNoiseRMS
+
+def Calc_Voltage_Noise_RMS(LSD, maxFreq):
+    """
+    See explanation of Calc_Power_Noise_RMS() and Lenny's Thesis 3.4.
+    """
+    PowerNoiseRMS = Calc_Power_Noise_RMS(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 / NumberofBits
+With 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 files
+f = 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 files
+data = 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/NoiseDensityHHCoilsFinalPowerSupplyWithPIBatteryInput.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
+
+samplingrateshielding, datashielding = wavfile.read("Z:/Users/Lenny/Power Supply/TimeSeriesOP549/FinalMeasurements/NoiseDensitySingleTestCoilWithPIBatteryInputWithShielding.wav")
+datashielding = datashielding*0.4/2**16
+datashieldingshaped = np.reshape(datashielding, [10,2000000])
+
+x_Final_av, y_Final_av = Average_Fourier_Uniform(dataFinalshaped, samplingrateFinal)
+x_shielding_av, y_shielding_av = Average_Fourier_Uniform(datashieldingshaped, samplingrateshielding)
+print(Calc_Voltage_Noise_RMS(x_Final_av, 25000))
+
+plt.rcParams["figure.figsize"] = [7,5]
+plt.rcParams["font.size"] = 12
+plt.figure(1, figsize=[12,7])
+plt.loglog(x_Final_av, y_Final_av, linewidth=1, label="Noise Final Powersupply in Carton Box", alpha=0.7)
+plt.loglog(x_shielding_av, y_shielding_av, linewidth=1, label="Noise Prototype with Shielding", 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("2 Measurements of 10M Samples at 6.52MHz averaged")
+plt.show()
+plt.close()
+
+