"""Analyze one uniformly sampled gyro axis; synthesize a clearly marked example.
python3 imu_noise.py --synthetic
python3 imu_noise.py stationary.csv
CSV columns: time_s,gyro_rad_s. Outputs go to the current directory.
"""
import argparse,json
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams.update({'svg.fonttype':'none','font.size':13})
def analyze(t,y):
    if len(t)<100 or not np.isfinite(t).all() or not np.isfinite(y).all():raise ValueError('Need at least 100 finite samples')
    d=np.diff(t);dt=float(np.median(d))
    if dt<=0 or np.any(d<=0) or np.max(np.abs(d-dt))>dt*.01:raise ValueError('Timestamps must increase with <=1% interval deviation; split gaps first')
    rows=[]
    for m in np.unique(np.logspace(0,np.log10(len(y)//10),35).astype(int)):
        k=len(y)//m;means=y[:k*m].reshape(k,m).mean(axis=1)
        rows.append((m*dt,np.sqrt(np.mean(np.diff(means)**2)/2),k-1))
    return {'samples':len(y),'dt_s':dt,'mean_rad_s':float(y.mean()),'std_rad_s':float(y.std(ddof=1))},np.array(rows)
if __name__=='__main__':
    p=argparse.ArgumentParser();p.add_argument('csv',nargs='?');p.add_argument('--synthetic',action='store_true');a=p.parse_args()
    if bool(a.csv)==a.synthetic:p.error('choose a CSV or --synthetic')
    if a.synthetic:
        rng=np.random.default_rng(7);t=np.arange(60000)/100;y=.01+rng.normal(0,.002,len(t));name='imu-synthetic'
        np.savetxt(name+'.csv',np.c_[t,y],delimiter=',',header='time_s,gyro_rad_s',comments='')
    else:
        data=np.genfromtxt(a.csv,delimiter=',',names=True);t=data['time_s'];y=data['gyro_rad_s'];name='imu-analysis'
    stats,rows=analyze(t,y);stats['synthetic']=a.synthetic;Path(name+'-stats.json').write_text(json.dumps(stats,indent=2)+'\n')
    np.savetxt(name+'-allan.csv',rows,delimiter=',',header='tau_s,allan_rad_s,adjacent_pairs',comments='')
    plt.figure(figsize=(9,5));plt.loglog(rows[:,0],rows[:,1],'o-',label='non-overlapping Allan deviation')
    if a.synthetic:plt.loglog(rows[:,0],.002/np.sqrt(rows[:,0]*100),'--',label='white-noise expectation')
    plt.xlabel('Averaging time tau (s)');plt.ylabel('Allan deviation (rad/s)');plt.title(('Synthetic: ' if a.synthetic else 'CSV input: ')+f'{len(y)} samples');plt.legend();plt.tight_layout()
    path=Path(name+'.svg');plt.savefig(path,metadata={'Date':None});s=path.read_text();i=s.index('<svg');j=s.index('>',i);s=s[:j]+' role="img" aria-labelledby="title desc"'+s[j:j+1]+f'<title id="title">{name} Allan deviation</title><desc id="desc">Non-overlapping means; units radians per second. Pair counts are in the CSV.</desc>'+s[j+1:];path.write_text(s)
    print(json.dumps(stats,indent=2))