Homec4science
Diffusion LC-Sampling-Theory (master) Browse

R12761/src/processing.py
0144bc3470bbmaster

Last Change
Raw File

processing.py
View Options

import numpy as np
import numpy.fft as fft
import jax.numpy as jnp
from jax import lax
INTERPOLATOR = 'sinc'
SIGNAL_TYPE = 'create' #'create': create a new emyulated ECG signal, 'load': load previous ECG signal, 'random': random sinusoidal mix signal
def get_max_freq(signal, freq, use_energy, stopping_prec):
dt = 1/freq
signal_0_mean = signal-np.mean(signal)
f_sig = fft.fft(signal_0_mean)[0:len(signal_0_mean)//2]
freq_sig = np.fft.fftfreq(signal.size, d=dt)[0:len(signal_0_mean)//2]
if use_energy:
threshold = sum(abs(f_sig)**2)/stopping_prec
eval_coeff = lambda x: abs(x)**2
else:
threshold = max(abs(f_sig))/stopping_prec
eval_coeff = lambda x: abs(x)
last_significant_coef_pos = 0
magn_f_c = 0
for i,f_c in enumerate(f_sig[-1::-1]):
magn_f_c += eval_coeff(f_c)
if magn_f_c > threshold:
last_significant_coef_pos = len(f_sig)-1-i
break
max_freq = freq_sig[last_significant_coef_pos]
return max_freq
def get_nyquist_freq_dataset(data,freq,use_energy,stopping_perc):
max_f = 0
for pt in data:
f = get_max_freq(pt,freq,use_energy,stopping_perc)
max_f = max(max_f,f)
return 2.1*max_f
def normalize(signal):
return signal/(jnp.max(signal)-jnp.min(signal))
def normalize_dataset(data):
out =[]
for pt in data:
out.append(normalize(pt))
return jnp.array(out)
# ----------------------------------
def gaussian(x, mu, sig):
return jnp.exp(-jnp.power(x - mu, 2.) / (2 * jnp.power(sig, 2.)))
def mix_gaussian_lvl_crossing(signal,mu_s,sigma):
return signal*jnp.sum(jnp.array([gaussian(signal,mu,sigma) for mu in mu_s]),axis = 0)
def lvl_crossing_stupid(x,mus):
out = jnp.zeros((len(x)-1,))
for i in range(1,len(x)):
for mu in mus:
#val = lax.cond((x[i-1]<mu and x[i]>=mu) or (x[i-1]>mu and x[i]<=mu),x[i],val)
if (x[i-1]<mu and x[i]>=mu) or (x[i-1]>mu and x[i]<=mu):
out = out.at[i].set(mu) #GRAD IS 0 IF I USE X[I] BUT NOT IF I USE MU
return out
def mix_gaussian_lvl_crossing_dataset(data,mu_s,sigma):
out = []
for signal in data:
out.append(mix_gaussian_lvl_crossing(signal,mu_s,sigma))
return np.array(out)
def sinc_interpolation_freq_parametrized(samples,freq,freq_in,time_base):
dt_samples = 1/freq_in
sinc_matr = jnp.array([jnp.sinc((time_base-n*dt_samples)*freq) for n in range(len(samples))])
return jnp.dot(samples,sinc_matr)
'''
out = samples[0]*jnp.sinc((time_base-0*dt_samples)*freq)
for n,f in enumerate(samples[1:]):
out += f*jnp.sinc((time_base-n*dt_samples)*freq) #interpolation with sync frequyency parametrized
return out
'''
def RMSE(v1,v2):
return jnp.sqrt(jnp.sum((v1-v2)**2)/len(v1))
# ----------------------------------
def sinc_interp(samples,freq_in,freq_out):
dt_samples = 1/freq_in
dt_new = 1/freq_out
new_time_base = np.arange(0,len(samples)*dt_samples,dt_new)
out = np.zeros((len(new_time_base)))
for n,f in enumerate(samples):
out += f*np.sinc((new_time_base-n*dt_samples)/dt_samples) #Whitaker Shannon
return out
def interpolate(signal,freq_in,freq_out,type = INTERPOLATOR):
if type == "sinc":
resampled = sinc_interp(signal,freq_in,freq_out)
else:
print("No interpolator recognized")
resampled = signal
return resampled
def interpolate_dataset(data,freq_in,freq_out,type = INTERPOLATOR):
out = []
for pt in data:
out.append(interpolate(pt,freq_in,freq_out))
return np.array(out)
c4science · Help