from __future__ import division
from deltasigma import *

OSR = 32
order = 5
H = synthesizeNTF(order, OSR, 1)

figure(figsize=(20, 4))
N = 8192
fB = int(np.ceil(N/(2.*OSR)))
ftest = np.floor(2./3.*fB)
u = 0.5*np.sin(2*np.pi*ftest/N*np.arange(N))
v, xn, xmax, y = simulateDSM(u, H)
t = np.arange(301)
step(t, u[t],'r')
hold(True)
step(t, v[t], 'g')
axis([0, 300, -1.2, 1.2])
xlabel('Sample Number')
ylabel('u, v')
title('Modulator Input & Output');

f = np.linspace(0, 0.5, N/2. + 1)
spec = np.fft.fft(v * ds_hann(N))/(N/4)
plot(f, dbv(spec[:N/2. + 1]),'b', label='Simulation')
figureMagic([0, 0.5], 0.05, None, [-120, 0], 20, None, (16, 6), 'Output Spectrum')
xlabel('Normalized Frequency')
ylabel('dBFS')
snr = calculateSNR(spec[2:fB+1], ftest - 2)
text(0.05, -10, 'SNR = %4.1fdB @ OSR = %d' % (snr, OSR), verticalalignment='center')
NBW = 1.5/N
Sqq = 4*evalTF(H, np.exp(2j*np.pi*f)) ** 2/3.
hold(True)
plot(f, dbp(Sqq * NBW), 'm', linewidth=2, label='Expected PSD')
text(0.49, -90, 'NBW = %4.1E x $f_s$' % NBW, horizontalalignment='right')
legend(loc=4);

snr_pred, amp_pred, _, _, _ = predictSNR(H, OSR)
snr, amp  = simulateSNR(H, OSR)

plot(amp_pred, snr_pred, '-', amp, snr, 'og-.')
figureMagic([-100, 0], 10, None, [0, 100], 10, None, (16, 6),'SQNR')
xlabel('Input Level (dBFS)')
ylabel('SQNR (dB)')
pk_snr, pk_amp = peakSNR(snr, amp)
text(-25, 85, 'peak SNR = %4.1fdB\n@ OSR = %d\n' % (pk_snr, OSR), horizontalalignment='right');

f0 = 1./8
OSR = 64
order = 8
H = synthesizeNTF(order, OSR, 1, 1.5, f0)

figure(figsize=(20, 4))
fB = int(np.ceil(N/(2. * OSR)))
ftest = int(np.round(f0*N + 1./3 * fB))
u = 0.5*np.sin(2*np.pi*ftest/N*np.arange(N))
v, xn, xmax, y = simulateDSM(u, H)
t = np.arange(301)
step(t, u[t], 'r')
hold(True)
step(t, v[t], 'g')
axis([0, 300, -1.2, 1.2])
xlabel('Sample Number')
ylabel('u, v')

spec = np.fft.fft(v*ds_hann(N))/(N/4)
plot(np.linspace(0, 0.5, N/2. + 1), dbv(spec[:N / 2 + 1]), label='Simulation')
figureMagic([0, 0.5], 0.05, None, [-140, 0], 20, None, (16, 6), 'Output Spectrum')
f1 = int(np.round((f0 - 0.25/OSR) * N))
f2 = int(np.round((f0 + 0.25/OSR) * N))
snr = calculateSNR(spec[f1:f2+1], ftest - f1)
text(0.15, -10, 'SNR = %4.1f dB @ OSR=%d)' % (snr, OSR), verticalalignment='center')
grid(True)
xlabel('Normalized Frequency')
ylabel('dBFS/NBW')
NBW = 1.5/N
Sqq = 4*evalTF(H, np.exp(2j*np.pi*f))**2/3.
hold(True)
plot(f, dbp(Sqq*NBW), 'm', linewidth=2, label='Expected PSD')
text(0.475, -90, 'NBW=%4.1E x $f_s$' % NBW, horizontalalignment='right', verticalalignment='center')
legend(loc=4);

snr_pred, amp_pred, _, _, _ = predictSNR(H, OSR, None, f0)
snr, amp = simulateSNR(H, OSR, None, f0)

plot(amp_pred, snr_pred, '-b', amp, snr, 'og-.')
figureMagic([-110, 0], 10, None, [0, 110], 10, None, (16, 6), 'SQNR')
xlabel('Input Level (dBFS)')
ylabel('SQNR (dB)')
pk_snr, pk_amp = peakSNR(snr, amp)
text(-20, 95, 'peak SNR = %4.1fdB\n@ OSR = %d\n' % (pk_snr, OSR), horizontalalignment='right');

figure(figsize=(16, 20))
colors = ('m', 'b')
Hinf_list = [2, 8]
i = -1
for col, Hinf in zip(colors, Hinf_list):
    i += 2
    OSR = 8
    M = 16
    H = synthesizeNTF(7, OSR, 1, Hinf)
    N = 8192
    fB = int(np.ceil(N/(2.*OSR)))
    ftest = int(np.floor(2./7*fB))
    u = 0.5*M*np.sin(2*np.pi*ftest/N*np.arange(N))
    v, xn, xmax, y = simulateDSM(u, H, M + 1)
    
    subplot(int('42' + str(i)))
    t = np.arange(101)
    step(t, u[t], 'g')
    hold(True)
    step(t, v[t], col)
    figureMagic([0, 100], 10, None, [-M, M], 2, None, None,'Input & Output')
    xlabel('Sample Number')
    ylabel('u, v')
    
    subplot(222)
    snr, amp = simulateSNR(H, OSR, None, 0., M + 1)
    plot(amp, snr,'o' + col, amp, snr, '--' + col)
    hold(True)
    figureMagic([-120, 0], 10, None, [0, 120], 10, None, None,'SQNR')
    xlabel('Input Level (dBFS)')
    ylabel('SNR (dB)')
    pk_snr, pk_amp = peakSNR(snr, amp)
    text(-13, pk_snr, 'peak SNR = %4.1fdB\n@ OSR = %d\n' % (pk_snr, OSR), horizontalalignment='right', color=col)
    
    subplot(212)
    f = np.linspace(0, 0.5, N/2. + 1)
    spec = np.fft.fft(v*ds_hann(N))/(M*N/4)
    plot(f, dbv(spec[:N/2 + 1]), col)
    snr = calculateSNR(spec[2:fB + 1], ftest - 2)
    text(0.1, 10*(i - 4), 'SNR = %4.1fdB @ OSR=%d' % (snr, OSR), verticalalignment='center', color=col)
    figureMagic([0, 0.5], 0.05, None, [-160, 0], 20, None, None,'Output Spectrum')
    xlabel('Normalized Frequency')
    ylabel('dBFS')
    
    NBW = 1.5/N
    Sqq = 4*evalTF(H, np.exp(2j*np.pi*f))**2/(3.*M**2)
    hold(True)
    plot(f, dbp(Sqq*NBW), 'c', linewidth=2)
    if i == 1:
        text(0.47, -110, 'NBW=%4.1E x $f_s$ '% NBW, horizontalalignment='right')
title('15-step 7th-order Lowpass')
tight_layout()