-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathtestScript.py
More file actions
197 lines (154 loc) · 5.39 KB
/
Copy pathtestScript.py
File metadata and controls
197 lines (154 loc) · 5.39 KB
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
import sys
import os
# Add build directory to sys.path
# Script is in examples/ directory, so need to go up one level to project root
script_dir = os.path.dirname(os.path.abspath(__file__))
project_root = os.path.dirname(script_dir)
build_dir = os.path.join(project_root, "build")
sys.path.insert(0, build_dir)
import wave
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import periodogram
from AmpClippers import AmpegB15Portaflex1A, SoldanoSuperLead100_2B
from ClippingCircuits import BiasDiodeGE_A
from EQCircuits import Baxandall3
from FilterCircuits import AllPassFilter
from Lite import ButterworthSallenKeyHPF_60Hz
from PedalCircuits import ProCoRat
from AmpTonestacks import FenderBassmanToneStack
# SET CONSTANTS
Fs = 48000 # Sample rate
Ts = 1/Fs # Period
figure, axis = plt.subplots(2, 2, tight_layout=True)
###########################
###### PLOT WAVEFORM ######
###########################
# Define Circuit
effect = AmpegB15Portaflex1A()
# effect = AllPassFilter()
# effect = SoldanoSuperLead100_2B()
# effect = BiasDiodeGE_A()
# effect = Baxandall3()
# effect = ButterworthSallenKeyHPF_60Hz()
# effect = ProCoRat()
# effect = FenderBassmanToneStack()
# Generate a waveform for testing
f = 500 # fundamental frequency
numCycles = 3
durSec = numCycles/f # duration of waveform
t = np.arange(0, durSec, Ts) # time vector
x1 = np.sin(2*np.pi*f*t).astype(np.float32)
# Set up an output variable
y1 = x1.astype(np.float32) # must cast to dtype=np.float32
numSamples = np.size(y1)
# Set parameters (if necessary)
newParams = effect.setParametersNoSmoothing([0.2, 0.5])
# Prepare the circuit with Fs and buffer size
effect.prepare(Fs, 1024)
# Process function
effect.process(x1, y1, numSamples, 0)
# Plot waveform
axis[0, 0].plot(t,x1,t,y1)
axis[0, 0].set_title('Waveform')
axis[0, 0].set_xlabel('Time (sec.)')
axis[0, 0].set_ylabel('Amplitude')
axis[0, 0].legend(['Input', 'Output'])
################################
###### FREQUENCY RESPONSE ######
################################
# Define Circuit
effect = AmpegB15Portaflex1A()
# Impulse response signal
x2 = scipy.signal.unit_impulse(4096).astype(np.float32)
numSamples = np.size(x2)
# Set parameters (if necessary)
newParams = effect.setParametersNoSmoothing([0.2, 0.5])
# Prepare the circuit with Fs and buffer size
effect.prepare(Fs, 1024)
# Process
y2 = x2.astype(np.float32)
effect.process(x2, y2, numSamples, 0)
# Compute frequency response
W, H = scipy.signal.freqz(y2, 1, worN=4096, fs=Fs)
Amp = 20*np.log10(abs(H) + 1e-10) # Convert to dB
# Plot frequency response
axis[0, 1].semilogx(W, Amp)
# plt.axvline(fc, color='k')
axis[0, 1].set_ylim(-20, max(Amp)+5)
axis[0, 1].set_xlim(20, 20000)
axis[0, 1].set_xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000], ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"])
axis[0, 1].set_title('Frequency Response')
axis[0, 1].set_xlabel('Freq. (Hz)')
axis[0, 1].set_ylabel('Amplitude (dB)')
#######################################
###### TOTAL HARMONIC DISTORTION ######
#######################################
# Define Circuit
effect = AmpegB15Portaflex1A()
# Generate input signal
f = 500 # fundamental frequency
numCycles = 200
durSec = numCycles/f
t = np.arange(0, durSec, Ts) # time vector
x3 = np.sin(2*np.pi*f*t).astype(np.float32) # sine wave
numSamples = np.size(x3)
# Set parameters (if necessary)
newParams = effect.setParametersNoSmoothing([0.2, 0.5])
# Prepare the circuit with Fs and buffer size
effect.prepare(Fs, 1024)
# Process
y3 = x3.astype(np.float32)
effect.process(x3, y3, numSamples, 0)
# Compute the periodogram
f, Pxx = periodogram(y3, Fs) # compute power spectral density
# Convert power spectral density (PSD) to dB
Pxx_db = 10 * np.log10(Pxx)
# Compute the total harmonic distortion (THD)
fundamental_idx = np.argmax(Pxx) # index of fundamental freq.
fundamental_power = Pxx[fundamental_idx] # power of fundamental freq.
harmonic_powers = Pxx[fundamental_idx*2::fundamental_idx] # power of harmonics
thd = np.sqrt(np.sum(harmonic_powers)) / fundamental_power # compute THD
# Downsample the periodogram for plotting
max_points = 10000 # Maximum number of data points for plotting
if len(f) > max_points:
factor = len(f) // max_points
f_downsampled = scipy.signal.resample(f, len(f) // factor)
Pxx_db_downsampled = scipy.signal.resample(Pxx_db, len(f) // factor)
else:
f_downsampled = f
Pxx_db_downsampled = Pxx_db
# Plot the periodogram
axis[1, 0].plot(f_downsampled, Pxx_db_downsampled)
axis[1, 0].set_xlabel('Frequency (Hz)')
axis[1, 0].set_ylabel('Power Spectral Density (dB)')
axis[1, 0].set_title('Periodogram with Harmonic Labels')
axis[1, 0].grid(True)
axis[1, 0].set_ylim(-75, 0)
# Print the THD value
print('Total Harmonic Distortion (THD):', thd)
######################
###### DC SWEEP ######
######################
# Define Circuit
effect = AmpegB15Portaflex1A()
minValue = -1
maxValue = 1
stepSize = 0.01
x = np.arange(minValue, maxValue, stepSize).astype(np.float32)
numSamples = np.size(x)
# Set parameters (if necessary)
newParams = effect.setParametersNoSmoothing([0.2, 0.5])
# Prepare the circuit with Fs and buffer size
effect.prepare(Fs, 1024)
# Process
y = x.astype(np.float32)
effect.processInPlace(y, numSamples, 0)
# Plot DC sweep
axis[1, 1].plot(x, x, x, y)
axis[1, 1].set_title('DC Sweep')
axis[1, 1].set_xlabel('Input Amplitude')
axis[1, 1].set_ylabel('Output Amplitude')
axis[1, 1].legend(['Input', 'Output'])
plt.show()