"""
This script contains the method *magic_function* which is used by
Y Feng et al. in their paper *General noise suppression scheme with
reference detection in heterodyne nonlinear spectroscopy*. It would be
interesting to implement this algorithm into the measurement software to
see if it can improve the signal to noise ratio etc.
The original matlab scripts from Feng as well as their paper were used
to translate the algorithm into Python.
"""
#%%
# Add parent directories to path for imports
import os, sys, inspect
currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
parentdir = os.path.dirname(currentdir)
sys.path.insert(0, parentdir)
sys.path.insert(1, os.path.join(parentdir, ""))
import numpy as np
import pprint as p
import h5py
from matplotlib import pyplot as plt
[docs]def magic_function(signal, reference):
# 1,32
avg_sig = np.average(signal.T, axis=1)
# print("avg_sig:", avg_sig.shape)
# print(avg_sig)
# 1, 32
avg_ref = np.average(reference.T, axis=1)
# print("avg_ref:", avg_ref.shape)
# Center data by element wise subtracting
# the average from the data
# shape (samples_to_acquire, no_pixels)
# 9922, 32
centered_sig = signal - avg_sig
centered_ref = reference - avg_ref
# Fit data
# Calculate covariance matrix of reference values
# This does not appear to be a real covariance #??????
# shape (no_pixels, no_pixels)
# 32, 32
cov_ref = np.matmul(centered_ref.T, centered_ref)
# using Cholesky factorization to inverse the
# covariance matrix
# shape (no_pixels, np_pixels)
inv_cov_ref = np.linalg.cholesky(cov_ref).T
# shape (no_pixels, no_pixels)
temp = np.matmul(centered_sig.T, centered_ref)
# shape (no_pixels, no_pixels)
# x*A=B is equivalent to A.'*x.'=B.'
# Solve the system of equations
Q = np.matmul(temp, np.linalg.inv(inv_cov_ref))
# Get the coefficient matrix
coeff = np.linalg.solve(inv_cov_ref, Q.T)
# Get the numbers of sample which were acquired
sample_acquired = centered_sig.shape[0]
# Calculate the variance of the signal
signal_variance = np.var(centered_sig, axis=0, ddof=1)
d = np.sum(Q ** 2, axis=1) / (sample_acquired - 1)
# Calculate standard deviation of the signal
# noise before referencing
signal_stabw = np.sqrt(signal_variance)
# noise after referencing
snr = np.sqrt(signal_variance - d.T)
return coeff, signal_stabw, snr
if __name__ == "__main__":
# Load in matlap file as arrays
# h5py is used because the .mat file
# is version 7.0
file = h5py.File("signal_reference_values.mat", "r")
file.keys()
# Create arrays dictionary
arrays = {}
for k, v in file.items():
arrays[k] = np.array(v)
# Read out reference and signal array
# from dictionary
# shape is (samples_to_acquire, no_pixels)
reference = arrays["ref"].T
signal = arrays["sig"].T
# print("reference:", reference.shape)
# Average over all laser shots
# shape (no_pixels,)
# avg_ref = np.average(reference.T, axis=1)
# avg_sig = np.average(signal.T, axis=1)
# Plot averaged signal and reference
# fig, axes = plt.subplots(nrows= 2, ncols=1)
# axes[0].plot(avg_ref.T)
# axes[1].plot(avg_sig.T)
# plt.show()
coeff, signal_stabw, snr = magic_function(signal, reference)