sinusoid recovery from low SNR signal

non_linear_filtering.py

import numpy as np
from scipy.signal import firwin, lfilter, butter, medfilt
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt


# sampling frequency in Hz
Fs = 8000

# center frequency of a band-limited signal
fc = 50

# Observe duration in seconds, forumulate a time-axis
T = 1
tvals = np.arange(0, T, 1/Fs)

fig=plt.figure(figsize=(10,10))
s_t = np.sin(2*np.pi*fc*tvals)

# Noise amplitude 5x signal strength
noiseA = 5
n_t = np.random.rand(len(s_t)) * noiseA * np.sqrt(np.power(10, -0.3))

plt.subplot(2,2,1)
plt.plot(tvals, s_t)
plt.ylabel('band-limited signal')
plt.grid()
# Signal + noise
x_t = s_t + n_t


plt.subplot(2,2,2)
plt.plot(tvals, x_t)
plt.ylabel('Low SNR noisy signal')
plt.grid()

# Simple FIR/IIR filtering, taps and filter order respectively
N = 31
L = 3
taps_lowpass = firwin(N, (fc)/Fs)
b,a = butter(L, 2 * (fc)/Fs, btype='low')
lp_filtered_signal = lfilter(taps_lowpass, 1.0, x_t)
bp_filtered_signal = lfilter(b, a, x_t)

plt.subplot(2,2,3)
plt.plot(tvals, lp_filtered_signal)
plt.grid()
plt.xlabel('Time')
plt.ylabel('FIR Low pass filtered signal')
plt.tight_layout()

plt.subplot(2,2,4)
plt.plot(tvals, bp_filtered_signal)
plt.xlabel('Time')
plt.ylabel('IIR Low pass filtered signal')
plt.tight_layout()
plt.grid()

plt.figure(2)
# Frequency analysis for non-linear filtering
X_w = np.fft.fft(x_t)
# Filter low energy noise components
max_freq_comp = np.max(np.abs(X_w))
X_w_ignore_low = np.where(np.abs(X_w) > max_freq_comp/10, X_w, 0+0j)
# Filter DC component
X_w_ignore_low[0] = 0+0j

F_w = np.fft.fftfreq(tvals.shape[-1])
y_t = np.fft.ifft(X_w_ignore_low)
plt.title('Non-linear filtering of AWGN tempered signal')
plt.subplot(2,2,1)
plt.plot(tvals, x_t)
plt.xlabel('Time')
plt.ylabel('Low SNR noisy signal')
plt.grid()

plt.subplot(2,2,2)
plt.plot(F_w, np.abs(X_w))
plt.xlabel('Frequency axis')
plt.ylabel('Frequency response noisy signal')
plt.grid()

plt.subplot(2,2,3)
plt.plot(F_w, np.abs(X_w_ignore_low))
plt.xlabel('Frequency axis')
plt.ylabel('Frequency response after non-linear filtering')
plt.grid()

plt.subplot(2,2,4)
plt.plot(tvals, y_t)
plt.ylabel('Recovered signal')
plt.xlabel('Time')
plt.grid()
plt.show()

Quick plotting python

import numpy as np
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt


fig = plt.figure(figsize=(10,10))
dimming_master_target=np.array([0, 1000, 2029, 2400, 2600, 2750, 3000, 3350, 3700, 4000],
				dtype=float)

# HMI coding 10%
tl_h010_b000=np.array([0.07534,0.07534,0.07534,0.2564,0.3588,0.4521,0.6628,0.9579,1,1])
sl_h010_b000=np.array([972, 972, 972, 1584, 1930, 2245, 2956, 3953, 4095, 4095])
tl_h010_b050=np.array([0.1753,0.1753,0.1753,0.3369,0.4281,0.5113,0.6993,0.9624,1,1])
sl_h010_b050=np.array([1310, 1310, 1310, 1856, 2164, 2445, 3079, 3968, 4095, 4095])
tl_h010_b100=np.array([0.2753,0.2753,0.2753,0.4173,0.4975,0.5706,0.7357, 0.967,1,1])
sl_h010_b100=np.array([1648, 1648, 1648, 2127, 2398, 2645, 3203, 3983, 4095, 4095])
plt.subplot(221)
plt.plot(dimming_master_target, tl_h010_b000, label="brightness=0")
plt.plot(dimming_master_target, tl_h010_b050, label="brightness=50")
plt.plot(dimming_master_target, tl_h010_b100, label="brightness=100")
plt.title("Target Luminance: HMI Coding 10%")
plt.xlabel('Dimming master target function light')
plt.grid()
plt.legend()

plt.subplot(222)
plt.plot(dimming_master_target, sl_h010_b000, label="brightness=0")
plt.plot(dimming_master_target, sl_h010_b050, label="brightness=50")
plt.plot(dimming_master_target, sl_h010_b100, label="brightness=100")
plt.title("Standardized Luminance: HMI Coding 10%")
plt.xlabel('Dimming master target function light')
plt.grid()
plt.legend()

# HMI coding 100%
tl_h100_b000=np.array([0,0, 0, 0, 0,0,0.3346,0.9168,1,1])
sl_h100_b000=np.array([718, 718, 718, 718, 718, 718, 1848, 3814, 4095, 4095])
tl_h100_b050=np.array([0.1753,0.1753,0.1753,0.3369,0.4281,0.5113,0.6993,0.9624,1,1])
sl_h100_b050=np.array([1310, 1310, 1310, 1856, 2164, 2445, 3079, 3968, 4095, 4095])
tl_h100_b100=np.array([1,1,1,1,1,1,1,1,1,1])
sl_h100_b100=np.array([4095, 4095, 4095, 4095, 4095, 4095, 4095, 4095, 4095, 4095])

plt.subplot(223)
plt.plot(dimming_master_target, tl_h100_b000, label="brightness=0")
plt.plot(dimming_master_target, tl_h100_b050, label="brightness=50")
plt.plot(dimming_master_target, tl_h100_b100, label="brightness=100")
plt.title("Target Luminance: HMI Coding 100%")
plt.xlabel('Dimming master target function light')
plt.grid()
plt.legend()

plt.subplot(224)
plt.plot(dimming_master_target, sl_h100_b000, label="brightness=0")
plt.plot(dimming_master_target, sl_h100_b050, label="brightness=50")
plt.plot(dimming_master_target, sl_h100_b100, label="brightness=100")
plt.title("Standardized Luminance: HMI Coding 100%")
plt.xlabel('Dimming master target function light')
plt.grid()
plt.legend()