-
Notifications
You must be signed in to change notification settings - Fork 0
/
second-order.py
69 lines (55 loc) · 1.48 KB
/
second-order.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
import ctcsound
import numpy as np
import matplotlib.pyplot as plt
def fft(x):
s = abs(np.fft.rfft(x)/(x.size/2))
s[0] *= 0.5 # DC has a factor of 2 in rfft()
s *= 1/max(s)
return 20*np.log10(s+0.000001)
cs = ctcsound.Csound()
cs.compile_(['','highorder.csd','-r5000000'])
cs.readScore('''
i1 0 1 1 500
i2 0 1 1 500
''')
sr = cs.sr()
chan = cs.nchnls()
buffer = np.zeros(int(sr)*chan)
r = cs.performKsmps()
cnt = 0
ksmps = cs.ksmps()
while r == 0:
if cnt < int(sr)*chan:
buffer[cnt:cnt+ksmps*chan] = cs.spout()
r = cs.performKsmps()
cnt += ksmps*chan
end = chan*3*round(sr/500)
freqs = np.arange(0,buffer.size/4+1)
time = np.arange(0,end//chan)/(sr)
fig,axs = plt.subplots(2,2)
axs[0,0].plot(time*1000, buffer[0:end:chan], 'k')
axs[0,1].plot(time*1000,buffer[1:end+1:chan], 'k')
axs[1,0].plot(freqs/1000, fft(buffer[0::chan]), 'k')
axs[1,1].plot(freqs/1000, fft(buffer[1::chan]), 'k')
axs[0,1].set_xlim(0,2.01)
axs[0,0].set_xlim(0,2.01)
axs[0,1].set_xlabel('time (ms)')
axs[0,0].set_xlabel('time (ms)')
axs[0,1].set_ylabel('amplitude')
axs[0,0].set_ylabel('amplitude')
axs[1,0].set_ylim(-90,0)
axs[1,1].set_ylim(-90,0)
axs[1,1].set_xlim(-1,20)
axs[1,0].set_xlim(-1,20)
axs[1,1].set_xlabel('frequency (KHz)')
axs[1,0].set_xlabel('frequency (KHz)')
axs[1,1].set_ylabel('magnitude (dB)')
axs[1,0].set_ylabel('magnitude (dB)')
axs[0,0].set_title('FM')
axs[0,1].set_title('PM')
for i in (0,1):
for j in (0,1):
axs[i,j].grid()
plt.tight_layout()
plt.savefig("fm-pm.png")
plt.show()