from __future__ import division
from __future__ import print_function

from math import ceil
from math import log
from math import sin
from math import pi
from fractions import gcd

import myhdl
print(myhdl.__version__)
from myhdl import *

def m_musicbox(clock, reset, note, nv, sample_rate=48e3, clock_rate=50e6):
    """ module to generate a "tone".
    
    Port Map
    --------
      clock : circuit synchronous clock
      reset : circult reset
      note  : digital signal for the note
      nv    : sample valid strobe
    """
    
    # Build the ROM table to hold the "note". 
    # Replace the following with which ever "note" 
    # sequence (algorithm) one desires
    Fs = sample_rate
    nmax = note.val.max
    f1,f2 = Fs/40, Fs/12
    print("The note has a tone at %.3f and %.3f" % (f1,f2,))
    nsmp = int((f1*f2)/gcd(f1,f2))
    note_rom = [(sin(p1)+sin(p2))/2 for p1,p2 in 
               zip([((2*pi)/Fs)*ii*f1 for ii in range(nsmp) ],
                   [((2*pi)/Fs)*ii*f2 for ii in range(nsmp)]) ]
    
    # convert the ROM from float to integer for the requested range
    note_rom = tuple([int(round(nmax*nn)) for nn in note_rom])
    
    # Signals and variable for the logic
    ticks_per_fs = int(ceil(clock_rate/sample_rate))
    sample_rate_cnt = Signal(intbv(1, min=0, max=ticks_per_fs+1))
    
    # note ROM current index
    Nsmp = len(note_rom)
    noteidx = Signal(intbv(0, min=0, max=Nsmp))
    
    @always_seq(clock.posedge, reset=reset)
    def rtl():        
        if sample_rate_cnt == ticks_per_fs:
            sample_rate_cnt.next = 1
            note.next = note_rom[noteidx]
            nv.next = True
            noteidx.next = noteidx + 1 if noteidx < Nsmp-1 else 0
        else:
            sample_rate_cnt.next = sample_rate_cnt + 1
            nv.next = False
        
    return rtl

# quick check
m_musicbox(Signal(bool(0)), ResetSignal(0, 0, False), 
           Signal(intbv(0, min=-16, max=16)), Signal(bool(0)) )

def test(Nsmp=10, Fs=48e3, Nmax=16, convert=False):
    
    Fclk = 4*48e3
    Ts = 1/Fs
    xv,tv = [],[]
    nmax = Nmax
    
    clock = Signal(bool(0))
    reset = ResetSignal(0, active=0, async=True)
    note = Signal(intbv(0, min=-nmax, max=nmax))
    nv = Signal(bool(0))
    
    # clock driver
    @always(delay(10))
    def tbclk():
        clock.next = not clock
        
    # reset driver
    def pulse_reset():
        reset.next = reset.active
        yield delay(100)
        reset.next = not reset.active
        yield clock.posedge
        
    # stimulus driver
    def _test():
        tbdut = m_musicbox(clock, reset, note, nv, Fs, Fclk)
                
        @instance
        def tbstim():
            t,scnt = 0,0
            yield pulse_reset()
            
            while scnt < Nsmp:
                if nv:
                    xv.append(int(note.val))
                    tv.append(t)
                    t += Ts
                    scnt += 1
                yield clock.posedge
                
            raise StopSimulation
            
        return tbclk, tbdut, tbstim
    
    # run the simulation, using _test as the stimulus
    Simulation(traceSignals(_test)).run()

    if convert:
        # convert the design to VHDL
        toVHDL(m_musicbox, clock, reset, note, nv,
               sample_rate=Fs, clock_rate=50e6)    
        
        toVerilog(m_musicbox, clock, reset, note, nv, 
                  sample_rate=Fs, clock_rate=50e6)    
    return tv,xv



# run the simulation
Fs = 48e3
Nmax = 2**(15-1)
tv,xv = test(Nsmp=100, Fs=Fs, Nmax=Nmax, convert=False)
print("Simulation complete")

# run the simulation and collect 1000 samples and plot
import matplotlib.pyplot as plt

tv,xv = test(Nsmp=Fs//200, Fs=Fs, Nmax=Nmax, convert=False)

fig,ax = plt.subplots(1, figsize=(9,3,))
ax.plot(tv, xv)
ax.set_ylim(-Nmax, Nmax)
ax.set_xlabel("Time [sec]")
ax.set_ylabel("Raw Amplitude")

# plot the power spectral density (see the spikes of the "tones")
fig,ax = plt.subplots(1, figsize=(9,3,))
tv,xv = test(Nsmp=Fs, Fs=Fs, Nmax=Nmax, convert=False)
p,f = ax.psd(xv, Fs=Fs, NFFT=8192)
ax.set_ylim(-60, 80)
ax.set_xlim(0, 6000)


# the following requires the myhdl_tools package 
# (http://www.bitbucket.org/cfelton/myhdl_tools)
from myhdl_tools import vcd

# rerun the test with less samples
test(Nsmp=5, Fs=Fs, Nmax=Nmax, convert=False)

# not a great VCD plotter but enough to get an idea,
# use gtkwave for real debug
vcd.parse_and_plot("_test.vcd")

# convert and run tools, this requires the gizflo
# package and can be found at http://www.github.com/cfelton/gizflo

tv,xv = test(Nsmp=2, Fs=Fs, Nmax=Nmax, convert=True)