1. 1
2. 2
3. 3
4. 4
5. 5
6. 6
7. 7
8. 8
9. 9
10. 10
11. 11
12. 12
13. 13
14. 14
15. 15
16. 16
17. 17
18. 18
19. 19
20. 20
21. 21
22. 22
23. 23
24. 24
25. 25
26. 26
27. 27
28. 28
29. 29
30. 30
31. 31
32. 32
33. 33
34. 34
35. 35
36. 36
37. 37
38. 38
39. 39
40. 40
41. 41
42. 42
43. 43
44. 44
45. 45
46. 46
47. 47
48. 48
49. 49
50. 50
51. 51
52. 52
53. 53
54. 54
55. 55
56. 56
57. 57
58. 58
59. 59
60. 60
61. 61
62. 62
63. 63
64. 64
65. 65
66. 66
67. 67
68. 68
69. 69
70. 70
71. 71
72. 72
73. 73
74. 74
75. 75
76. 76
77. 77
78. 78
79. 79
80. 80
81. 81
82. 82
83. 83
84. 84
85. 85
86. 86
87. 87
88. 88
89. 89
90. 90
91. 91
92. 92
93. 93
94. 94
95. 95
96. 96
97. 97
98. 98
99. 99
100. 100
101. 101
102. 102
103. 103
104. 104
105. 105
106. 106
107. 107
108. 108
109. 109
110. 110
111. 111
112. 112
113. 113
114. 114
115. 115
116. 116
117. 117
118. 118
119. 119
120. 120
121. 121
122. 122
123. 123
124. 124
125. 125
126. 126
127. 127
128. 128
129. 129
130. 130
131. 131
132. 132
133. 133
134. 134
135. 135
136. 136
137. 137

import bisect
import math
import struct
import sys

# Number of times to sample each second
bitrate = 44100

intro = [
    (1, 3, 3),
    (1, 3, 7),
    (1, 3, 10),
    (3, 4, 2),
    (1, 3, 1),
    (1, 3, 5),
    (1, 3, 8),
    (3, 4, 0),
    (1, 2, 11),
    (1, 3, 3),
    (1, 3, 6),
    (2, 3, 10),
    (1 / 2, 3, 9, 1 / 2),
    (1 / 2, 3, 8, 1 / 2),
    (6, 3, 7, 1 / 2),
]

melody = [
    (3, 3, 7),
    (3, 3, 3),
    (3, 3, 8),
    (2, 3, 7),
    (2, 3, 5),
    (3, 3, 7),
    (2, 3, 3),
    (2, 2, 10),
    (2, 3, 0),
    (6, 3, 3),
]

bass = [
    (1, 0, 7),
    (1, 2, 7),
    (1, 2, 7),
    (1, 0, 7),
    (1, 2, 3),
    (1, 2, 3),
    (1, 0, 7),
    (1, 2, 8),
    (1, 2, 8),
    (1, 0, 7),
    (1, 2, 7),
    (1, 0, 7),
    (1, 2, 5),
    (1, 0, 7),
    (1, 2, 7),
    (1, 2, 7),
    (1, 0, 7),
    (1, 2, 3),
    (1, 0, 7),
    (1, 1, 10),
    (1, 0, 7),
    (1, 2, 0),
    (1, 0, 7),
    (1, 2, 3),
    (1, 2, 3),
    (1, 0, 7),
    (1, 2, 3),
    (1, 2, 3),
]

music = []


def process(notes, start, speed=1, gain=1):
    """
    Adds a list of notes to the music list
    """
    t = start
    for note in notes:
        vol = 1
        if len(note) == 4:
            vol = note[3]
        music.append((t, t + note[0] / speed, note[1], note[2], vol * gain))
        t += note[0] / speed


# Process all lists of notes
process(intro, 0, 3)
process(melody, 8, 3)
process(bass, 8, 3, 1 / 2)
# And sort them
music.sort()


def freq(octave, step):
    """
    Returns the frequency of a note
    """
    return 55 * 2 ** (octave + step / 12 - 1)


def tone(f, t):
    """
    Returns the intensity of a tone of frequency f sampled at time t
    """
    # https://dsp.stackexchange.com/questions/46598/mathematical-equation-for-the-sound-wave-that-a-piano-makes
    # https://youtu.be/ogFAHvYatWs?t=254
    # return int(2**13*(1+square(t, 440*2**(math.floor(5*t)/12))))
    # Y = sum([math.sin(2 * i * math.pi * t * f) * math.exp(-0.0004 * 2 * math.pi * t * f) / 2**i for i in range(1, 4)])
    # Y += Y * Y * Y
    # Y *= 1 + 16 * t * math.exp(-6 * t)
    w = 2 * math.pi * f
    Y = 0.6 * math.sin(w * t) * math.exp(-0.0015 * w * t)
    Y += 0.4 * math.sin(2 * w * t) * math.exp(-0.0015 * w * t)
    Y += Y * Y * Y
    Y *= 1 + 16 * t * math.exp(-6 * t)
    return Y


def at(t):
    """
    Returns the total intensity of music sampled at time t
    """
    i = bisect.bisect(music, (t, 10**9))
    # This is actually pretty efficient ngl
    # Because people usually don't have that many overlapping notes
    ret = 0
    for j in range(max(i - 20, 0), i):
        m = music[j]
        if m[0] + m[1] >= t:
            ret += m[4] * tone(freq(m[2], m[3]), t - m[0])
    return int(2**12 * ret)


# Print out music encoded in s16 to standard output
for i in range(0 * bitrate, 16 * bitrate):
    sys.stdout.buffer.write(struct.pack("h", at(i / bitrate)))