-
Notifications
You must be signed in to change notification settings - Fork 1
/
dds_figure5.py
141 lines (108 loc) · 3.75 KB
/
dds_figure5.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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 21 14:15:39 2020
@author: DM1CR
"""
# python port of illustration skripts for
# Circuit cellar #217 - The darker side : DDS generator
# Robert Lacoste
import numpy as np
import matplotlib.pyplot as plt
# Global parameters
phaseregbitcount=16
clockfreq=1e6
dacbitcount=8
simulationsteps=1000
nyquistbandsforspectrumcalculation=6
#-----------------------------
# First simulation
phaseincrement=5169
# Deducted parameters
phaseregmax=2**phaseregbitcount
dacmax=2**(dacbitcount-1)
timeincrement=1/clockfreq
# steps vector
steps=np.arange(0,simulationsteps)
# Time vector
time=steps*timeincrement
# Phase vector
phase=np.mod(steps*phaseincrement,phaseregmax)
# Output sine lookup vector
scaledphase=phase*2*np.pi/phaseregmax
fout=(dacmax*np.sin(scaledphase)).astype(int)
# sine vector with better time resolution
fout_ext=np.arange(0,simulationsteps*nyquistbandsforspectrumcalculation)
for i in range(0,simulationsteps*nyquistbandsforspectrumcalculation-1):
fout_ext[i]=fout[int((i-1)/nyquistbandsforspectrumcalculation)]
# calculate spectrum
spectrum=np.fft.fft(fout_ext)
spectrum=2*spectrum/(simulationsteps*nyquistbandsforspectrumcalculation)
# However as the input signal is real then only the first half of the FFT
# is actually useful (the last half is the complex conjugate of the first half) :
usefulspectrum=spectrum[0:int(len(spectrum)/2)]
spectrumamplitudes=np.abs(usefulspectrum)
spectrumamplitudes[0]=dacmax #just to scale the graph
# sin(x)/x reference
x=np.arange(0, int(nyquistbandsforspectrumcalculation*simulationsteps/2))
x=x*np.pi/simulationsteps
sinxoverx=x
for i in range(1,len(x)-1):
sinxoverx[i]=dacmax*abs( np.sin(x[i])/x[i])
# Plot signals
plt.subplot(3,2,1)
plt.plot(phase[0:49])
plt.title('phase (N=16, P=16, B=8, W=5169)')
plt.subplot(3,2,3)
plt.plot(fout_ext[0:50*nyquistbandsforspectrumcalculation-1])
plt.title('Fout')
plt.subplot(3,2,5)
plt.plot(np.arange(0,len(spectrumamplitudes)), spectrumamplitudes,'k')
plt.plot(np.arange(0,len(spectrumamplitudes)), sinxoverx,'g')
plt.title('spectrum')
#-----------------------------
# second simulation
phaseincrement=15673
# Deducted parameters
phaseregmax=2**phaseregbitcount
dacmax=2**(dacbitcount-1)
timeincrement=1/clockfreq
# steps vector
steps=np.arange(0,simulationsteps)
# Time vector
time=steps*timeincrement
# Phase vector
phase=np.mod(steps*phaseincrement,phaseregmax)
# Output sine lookup vector
scaledphase=phase*2*np.pi/phaseregmax
fout=(dacmax*np.sin(scaledphase)).astype(int)
# sine vector with better time resolution
fout_ext=np.arange(0,simulationsteps*nyquistbandsforspectrumcalculation)
for i in range(0, simulationsteps*nyquistbandsforspectrumcalculation-1):
fout_ext[i]=fout[int((i-1)/nyquistbandsforspectrumcalculation)]
# calculate spectrum
spectrum=np.fft.fft(fout_ext)
spectrum=2*spectrum/(simulationsteps*nyquistbandsforspectrumcalculation)
# However as the input signal is real then only the first half of the FFT
# is actually useful (the last half is the complex conjugate of the first half) :
usefulspectrum=spectrum[0:int(len(spectrum)/2)]
spectrumamplitudes=np.abs(usefulspectrum)
spectrumamplitudes[0]=dacmax #just to scale the graph
# sin(x)/x reference
x=np.arange(0,int(nyquistbandsforspectrumcalculation*simulationsteps/2))
x=x*np.pi/simulationsteps
sinxoverx=x
for i in range(1,len(x)-1):
sinxoverx[i]=dacmax*abs(np.sin(x[i])/x[i])
# Plot signals
plt.subplot(3,2,2)
plt.plot(phase[0:49])
plt.title('phase (N=16, P=16, B=8, W=15673)')
plt.subplot(3,2,4)
plt.plot(fout_ext[0:50*nyquistbandsforspectrumcalculation-1])
plt.title('Fout')
plt.subplot(3,2,6)
plt.plot(np.arange(0,len(spectrumamplitudes)), \
spectrumamplitudes,'k', \
np.arange(0,len(spectrumamplitudes)), \
sinxoverx,'g')
plt.title('spectrum')