import numpy as np

# Volledig synthetische spanning/optische amplitude, geen opname van een mens.
fs=250.; t=np.arange(0.,10.,1/fs)
R=np.array([.50,1.30,2.12,2.90,3.72,4.52,5.31,6.14,6.93,7.75,8.55,9.34])
vertraging=.22+.020*np.sin(np.arange(len(R))*.9)
P=R+vertraging
ecg=np.zeros_like(t); ppg=np.zeros_like(t)
def gauss(mu,s): return np.exp(-.5*((t-mu)/s)**2)
for r,p in zip(R,P):
    ecg+=.12*gauss(r-.18,.035)-.15*gauss(r-.035,.012)
    ecg+=gauss(r,.012)-.25*gauss(r+.03,.015)+.30*gauss(r+.24,.065)
    ppg+=gauss(p,.055)+.25*gauss(p+.16,.045)
artefact=ppg+.95*gauss(5.12,.025)
def pieken(x,drempel,afstand=.35):
    # Drempels uitsluitend voor deze ontworpen golfvorm; geen medische detector.
    kandidaten=np.flatnonzero((x[1:-1]>x[:-2])&(x[1:-1]>=x[2:])&(x[1:-1]>drempel))+1
    over=[]
    for j in kandidaten:
        if not over or t[j]-t[over[-1]]>=afstand: over.append(j)
    return np.array(over,int)
ir=pieken(ecg,.60); ip=pieken(ppg,.60); ia=pieken(artefact,.60)
rr=np.diff(t[ir]); pp=np.diff(t[ip])
def cijfers(dt):
    return 60/np.mean(dt),1000*np.sqrt(np.mean(np.diff(dt)**2))
print("synthetische piekaantallen ECG/PPG/artefact",len(ir),len(ip),len(ia))
print("ECG hartfrequentie/RMSSD(ms)",np.round(cijfers(rr),6))
print("PPG pulsfrequentie/RMSSD(ms)",np.round(cijfers(pp),6))
print("artefact gemiddelde pulsfrequentie",cijfers(np.diff(t[ia]))[0])
# Onafhankelijke synthetische bewegingsvlag; niet uit de pulsfrequentie verzonnen.
kwaliteit=~((t>=5.00)&(t<=5.20))
goed=ia[kwaliteit[ia]]
print("na bewegingsvlag pieken",len(goed),"frequentie",cijfers(np.diff(t[goed]))[0])
print("perfecte PP-RR verschil max ms",1000*np.max(np.abs(np.diff(P)-np.diff(R))))
assert len(ir)==12 and len(ip)==12 and len(ia)==13 and len(goed)==12
assert np.allclose(np.diff(P),np.diff(R)+np.diff(vertraging))