"""Synthetic harmonic/filter study. No speech synthesis or clinical predictions.
Reproduce: python reproduce.py (requires matplotlib).
"""
import csv,json,math
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
ROOT=Path(__file__).resolve().parent

def transfer(f,fc):return 1/math.sqrt(1+(f/fc)**2)
def main():
    rows=[]
    for f0 in (100,200):
        for fc in (500,1500):
            for gain in (0,6):
                for n in range(1,21):
                    f=n*f0;a=1/n;h=transfer(f,fc);g=10**(gain/20)
                    rows.append({'fundamental_hz':f0,'harmonic_number':n,'frequency_hz':f,'cutoff_hz':fc,'gain_db':gain,'input_amplitude_au':a,'filter_magnitude':h,'filter_phase_rad':-math.atan(f/fc),'output_amplitude_au':a*h*g})
    single=[]
    for gain in (0,6):
        for n in range(1,21):
            a=1 if n==1 else 0
            single.append({'fundamental_hz':100,'harmonic_number':n,'frequency_hz':n*100,'cutoff_hz':500,'gain_db':gain,'input_amplitude_au':a,'output_amplitude_au':a*transfer(n*100,500)*10**(gain/20)})
    for filename,data in [('harmonics.csv',rows),('single-tone.csv',single)]:
        with (ROOT/filename).open('w',newline='',encoding='utf-8') as out:
            w=csv.DictWriter(out,fieldnames=list(data[0]));w.writeheader();w.writerows(data)
    checks={'gain_6db_amplitude_ratio':10**(6/20),'gain_ratio_max_error':0,'cutoff_attenuation_db':20*math.log10(transfer(500,500)),'cutoff_expected_db':-10*math.log10(2),'absent_components_remain_zero':all(q['output_amplitude_au']==0 for q in single if q['harmonic_number']>1),'frequency_retention':all(q['frequency_hz']==q['fundamental_hz']*q['harmonic_number'] for q in rows)}
    errors=[]
    for q in rows:
        if q['gain_db']==6:
            zero=next(z for z in rows if z['gain_db']==0 and z['fundamental_hz']==q['fundamental_hz'] and z['cutoff_hz']==q['cutoff_hz'] and z['harmonic_number']==q['harmonic_number'])
            errors.append(abs(q['output_amplitude_au']/zero['output_amplitude_au']-10**(.3)))
    checks['gain_ratio_max_error']=max(errors)
    assert checks['gain_ratio_max_error']<1e-12
    assert abs(checks['cutoff_attenuation_db']-checks['cutoff_expected_db'])<1e-12
    assert checks['absent_components_remain_zero'] and checks['frequency_retention']
    info={'title':'Synthetic harmonic source through a first-order low-pass filter','status':'Analytical illustrative study, not recorded speech, intelligibility or clinical evidence','parameters':{'fundamentals_hz':[100,200],'harmonics':list(range(1,21)),'source_amplitude':'1/n arbitrary units','cutoffs_hz':[500,1500],'gains_db':[0,6]},'equations':{'source':'s(t)=sum(n=1..20) sin(2*pi*n*f0*t)/n','filter':'H(f)=1/(1+j*f/fc)','magnitude':'1/sqrt(1+(f/fc)^2)','phase_rad':'-atan(f/fc)','output_amplitude':'(1/n)*abs(H(n*f0))*10^(gain_db/20)'},'rows':{'harmonics':len(rows),'single_tone':len(single)},'checks':checks,'limitations':['Arbitrary synthetic coefficients; no measured voice or vocal-tract response.','First-order electrical low-pass example is not a vocal-tract or electrolarynx model.','Periodic harmonics omit unvoiced speech and time-varying articulation.','A linear time-invariant filter and gain scale existing components; no missing frequency components are recovered.','No word-recognition, naturalness or clinical outcome is inferred.'],'source':{'url':'https://www.nidcd.nih.gov/news/multimedia/how-does-human-body-produce-voice-and-speech-text-version','verified':'2026-09-27','scope':'Official explanation of airflow, vocal-fold vibration and articulation; it does not specify or validate these model coefficients.'}}
    (ROOT/'study.json').write_text(json.dumps(info,indent=2)+'\n',encoding='utf-8')
    plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'axes.grid':True,'grid.alpha':.2})
    fig,axs=plt.subplots(1,3,figsize=(13,4.3),layout='constrained')
    for ax,f0,panel in zip(axs[:2],(100,200),('A','B')):
        ax.vlines([n*f0 for n in range(1,21)],0,[1/n for n in range(1,21)],color='#aaaaaa',alpha=.5,lw=4,label='Input: amplitude 1/n')
        for fc,color,marker in ((500,'#285ec7','o'),(1500,'#ad367f','s')):
            r=[q for q in rows if q['fundamental_hz']==f0 and q['cutoff_hz']==fc and q['gain_db']==0]
            ax.plot([q['frequency_hz'] for q in r],[q['output_amplitude_au'] for q in r],marker=marker,ms=3,lw=.9,color=color,label=f'Cutoff {fc} Hz, 0 dB gain')
        ax.set(title=f'{panel}. Fundamental {f0} Hz',xlabel='Frequency (Hz)',ylabel='Component amplitude (arbitrary units)',ylim=(0,1.05));ax.legend(fontsize=7)
    ax=axs[2]
    for gain,color,offset in ((0,'#285ec7',-10),(6,'#ad367f',10)):
        r=[q for q in single if q['gain_db']==gain]
        ax.vlines([q['frequency_hz']+offset for q in r],0,[q['output_amplitude_au'] for q in r],color=color,lw=2,label=f'{gain} dB gain')
        ax.plot([q['frequency_hz']+offset for q in r],[q['output_amplitude_au'] for q in r],'.',color=color)
    ax.set(title='C. Gain cannot restore absent components',xlabel='Frequency (Hz)',ylabel='Component amplitude (arbitrary units)',ylim=(-.05,2.1));ax.legend(fontsize=8)
    ax.text(.98,.55,'Single 100 Hz input tone\n500 Hz filter cutoff\nMarkers offset ±10 Hz for clarity',transform=ax.transAxes,ha='right',fontsize=8)
    fig.suptitle('Synthetic source and linear filtering • no speech or clinical outcomes',fontsize=12)
    for ext in ('png','svg','pdf'):fig.savefig(ROOT/f'voice-study.{ext}',dpi=200,bbox_inches='tight')
    plt.close(fig)
    print(json.dumps(checks,indent=2))
if __name__=='__main__':main()