Python assignment2 at University of Birmingham

ASSIGNMENT 2:System Identification

In this assignment you are presented with a number of ‘black box’ components, each of which contain code that models a certain electronic system. By investigating the response of each box to generated input (voltage) signals, you will be able to determine the type of system each box contain and estimate its parameters (e.g. resistance and inductance).
Before you continue with this assignment, please read Yr2_Computing_System_Identification.pdf and work though the example notebook black_box_example.ipynb first.

As before, you must store your student ID in the variable studentID and execute the cell below before starting the assignment. This will import the necessary Boxes objects. As in assignment 1, the final cell contains a function call (Boxes.check()) which you should run to check that the naming of variables and functions in each task is correct.

1
2
3
4
5
6
7
8
9
10
import numpy as np
import matplotlib.pyplot as plt

# Do not alter any of the code within this cell other than
# the value of studentID
from module_engine.assignment import Boxes
studentID = # Enter your student ID here

# Creating Box objects to be investigated in this assignment
t1_box1, t1_box2, t1_box3, t2_box, t3_box = Boxes.get_boxes(studentID)

Task 1

Marks available: 4

Each of the following circuits is contained in one of the boxes t1_box1, t1_box2 and t1_box3:

  • RL high pass filter
  • RC low pass filter
  • RLC band pass filter

By analysing the action of each box on custom signals and employing the .fft and .fftfreq functions from numpy.fft, identify the type of circuit contained in each box. A blank cell has been provided for your working, but you are free to insert as many additional cells as you require.

After completing your investigation set the variables RL_circuit, RC_circuit and RLC_circuit in the next cell to a value of 1, 2 or 3 to indicate which box matches the relevant circuit. For example, write RC_circuit = 1 if you believe t1_box1 to be the RC low pass filter or
RC_circuit = 2 if you believe that system is in t1_box2.

Reminder:
Use the .process() method to interact with the Boxes objects, as demonstrated in black_box_example.ipynb

Hint:
In order to obtain an informative result, your time domain input signal should contain frequency components ranging from 10 to 4000 Hz.

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
# Working for Task 1 (you are free to insert additional
# cells above/below this one if desired)
%matplotlib inline

def swept_sine(f1, f2, t):
"""
create a sine wave with a frequency sweeping between f1 and f2
f1: start frequency [Hz]
f2: stop frequency [Hz]
t: np.array of linearly spaced points in time [s]
"""
f = np.linspace(f1,f2,t.size)
output = np.sin(2*np.pi*f*t)
return output

# Sampling a signal correctly is a bit tricky, you need
# to make sure you know the total time, timestep and number
# of points accuractely. If in doubt, compute these from the
# time array you created and print the result.

def get_time_array(N=1000, duration=10, endpoint=False):
# N: Number of sample points
# duration: end point of the open interval

t = np.linspace(0, duration, num=N, endpoint=endpoint)
total_time = t[-1]-t[0] # different from duration!!
timestep = total_time/(N-1)
print("Data points = {}".format(N))
print("Total time = {} s".format(total_time))
print("Timestep = {} s".format(timestep))
print("Sampling Frequency = {} Hz".format((1/timestep)))

return t

# generate time array
t = get_time_array(N=40000, duration = 10.0)

# generate array containing input signal (voltage)
s_in = swept_sine(10,4000,t)

Data points = 40000
Total time = 9.99975 s
Timestep = 0.00025 s
Sampling Frequency = 4000.0 Hz

1
2
3
4
5
6
7
8
9
fig = plt.figure(figsize=(15,3))

# Plot input voltage as a function of time
plt.plot(t, s_in, 'b-')
plt.title('Input Signal')
plt.xlabel(r'$t$')
plt.ylabel(r'$V_{in}$')
plt.axis([0, 5, np.amin(s_in), np.amax(s_in)])
plt.show()

mark

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
# Pass input signal through the example_box
# and create an output signal (voltage):
s_out1 = t1_box1.process(t, s_in)
s_out2 = t1_box2.process(t, s_in)
s_out3 = t1_box3.process(t, s_in)

# And plot the result, comparing input to output
fig = plt.figure(figsize=(15,8))
# Plot input voltage as a function of time
plt.subplot(411)
plt.plot(t, s_in, 'b-')
plt.title('Input Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{in}$', fontsize=14)
plt.axis([0, 5, np.amin(s_in), np.amax(s_in)])

# Plot output voltage as a function of time
plt.subplot(412)
plt.plot(t, s_out1, 'r-')
plt.title('Output Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{out1}$', fontsize=14)
plt.axis([0, 5, np.amin(s_out1), np.amax(s_out1)])

# Plot output voltage as a function of time
plt.subplot(413)
plt.plot(t, s_out2, 'g-')
plt.title('Output Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{out2}$', fontsize=14)
plt.axis([0, 5, np.amin(s_out2), np.amax(s_out2)])

# Plot output voltage as a function of time
plt.subplot(414)
plt.plot(t, s_out3, 'y-')
plt.title('Output Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{out3}$', fontsize=14)
plt.axis([0, 5, np.amin(s_out3), np.amax(s_out3)])

# Automatically adjust subplots to fit figure area
plt.tight_layout()
# Display plots
plt.show()

mark

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
# Perform the Fourier Tranform on both signals
S_in = np.fft.fft(s_in)
S_out1 = np.fft.fft(s_out1)
S_out2 = np.fft.fft(s_out2)
S_out3 = np.fft.fft(s_out3)

# Next we generate sample points in the frequency domain,
# for this we need to know the duration and time step
# of the time array.
def get_frequency_array(t):
N = len(t)
timestep = np.diff(t)[0]
f = np.fft.fftfreq(N, d=timestep)
return f

f = get_frequency_array(t)

# We only want to plot the positive frequencies from the FFT
# this means we can ignore half of each array.
# Check the numpy.fftfreq documentation to see that
# frequencies are arranged as follows:
# f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even
# so we are interested in the elements f[1:n/2].
def get_positive_frequencies(A):
N=len(A)
if len(A)//2==0:
N2=int(N/2)
else:
N2=int((N+1)/2)
return A[1:N2]

f = get_positive_frequencies(f)
S_in = get_positive_frequencies(S_in)
S_out1 = get_positive_frequencies(S_out1)
S_out2 = get_positive_frequencies(S_out2)
S_out3 = get_positive_frequencies(S_out3)

# Creating some custom subplots
fg, ax = plt.subplots(4,1, figsize=(15,8), sharex=True)

# Plot the absolute values of transformed_output
ax[0].plot(f, np.abs(S_in), 'b-', label='input signal')
ax[0].plot(f, np.abs(S_out1), 'r-', label='output1 signal')
ax[0].plot(f, np.abs(S_out2), 'g-', label='output2 signal')
ax[0].plot(f, np.abs(S_out3), 'y-', label='output3 signal')
ax[0].legend()

ax[1].plot(f, np.abs(S_out1/S_in), 'k:', label='transfer function (output/input)')
ax[1].set_ylabel('|T| [V/V]')
ax[1].legend()

ax[2].plot(f, np.abs(S_out2/S_in), 'k:', label='transfer function (output/input)')
ax[2].set_ylabel('|T| [V/V]')
ax[2].legend()

ax[3].plot(f, np.abs(S_out3/S_in), 'k:', label='transfer function (output/input)')
ax[3].set_xlabel('Frequency [Hz]')
ax[3].set_ylabel('|T| [V/V]')
ax[3].legend()

ax[0].set_title('FFT of in- and output signals')
ax[0].set_ylabel('|Amplitude| [V/Hz]')

# Fine-tune figure; make subplots close to each other
# and hide x ticks for top plot.
fg.subplots_adjust(hspace=0)
plt.setp(ax[0].get_xticklabels(), visible=False)
plt.show()

mark

1
2
3
4
5
6
7
# Now, assign each of these variables a (unique) value 
# of 1, 2 or 3. For example if you know that t1_box2
# contains the RC_circuit, set RC_circuit=2, if instead
# you found the RC_circuit in Box t1_box1, set RC_circuit=1.
RL_circuit = 3
RC_circuit = 1
RLC_circuit = 2

Task 2

Marks available: 3

You are to work with t2_box in this task (defined in the first cell of the notebook). It is known that this box contains an RL high pass filter circuit, but with parameter values different from those of the RL high pass circuit in Task 1.

Employing a similar method to that used in Task 1, examine the characteristic transfer function of this circuit to determine an approximate value for the associated corner frequency, $f_0$. This corresponds to the frequency $f_0$ such that $|T(f_0)| =\sqrt{1/2}$.

The examples in Yr2_Computing_System_Identification.pdf refer to $T(\omega)$ while we sometimes use $T(f)$. These are the same function just plotted over a different axis ($\omega = 2\pi f$). It is easy to make mistakes when applying this scaling. You should always check if you are using the right frequency scaling for a given task.

After you complete you investigation store the value for f_0 in the variable corner_frequency. During marking, your value will be compared to the true value with a tolerance of 10 Hz. Ensure that your sample density in the frequency domain is sufficient to achieve this level of accuracy.

Furthermore, assuming that the resistance R of the circuit is 200 Ohms, calculate the inductance of the inductor in Henrys and store this in the variable inductance. Do not round this value.

Hints:

  • Your time domain input signal should contain frequency components ranging from 10 to 3000 Hz in order to obtain an informative result.
  • For finding a specific value in a numpy array you may want to use the numpy.argmin function.
    Play with and understand the following code:
1
2
3
4
a=np.array([10, 11, 12, 15, 20, 30])
index=np.argmin(np.abs(a-15))
print(index)
print(a[index])
  • To visually check whether you found the corner frequency correctly you can plot additional lines into the graph with the transfer fucntion, for example using plt.axvline(x=f0) and plt.axhline(y=1/np.sqrt(2),linestyle='--').
1
2
3
4
5
6
7
# Define parameters for Task 2 - explore frequencies up to 3000Hz

# generate time array
t = get_time_array(N=30000, duration = 10.0)

# generate array containing input signal (voltage)
s_in = swept_sine(10,3000,t)

Data points = 30000
Total time = 9.999666666666666 s
Timestep = 0.0003333333333333333 s
Sampling Frequency = 3000.0 Hz

1
2
3
4
5
6
7
8
9
fig = plt.figure(figsize=(15,3))

# Plot input voltage as a function of time
plt.plot(t, s_in, 'b-')
plt.title('Input Signal')
plt.xlabel(r'$t$')
plt.ylabel(r'$V_{in}$')
plt.axis([0, 5, np.amin(s_in), np.amax(s_in)])
plt.show()

mark

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
# Pass input signal through the example_box
# and create an output signal (voltage):
s_out = t2_box.process(t, s_in)

# And plot the result, comparing input to output
fig = plt.figure(figsize=(15,8))
# Plot input voltage as a function of time
plt.subplot(411)
plt.plot(t, s_in, 'b-')
plt.title('Input Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{in}$', fontsize=14)
plt.axis([0, 5, np.amin(s_in), np.amax(s_in)])

# Plot output voltage as a function of time
plt.subplot(412)
plt.plot(t, s_out, 'r-')
plt.title('Output Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{out}$', fontsize=14)
plt.axis([0, 5, np.amin(s_out), np.amax(s_out)])

# Automatically adjust subplots to fit figure area
plt.tight_layout()
# Display plots
plt.show()

mark

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
# Perform the Fourier Tranform on both signals
S_in = np.fft.fft(s_in)
S_out = np.fft.fft(s_out)

f = get_frequency_array(t)

f = get_positive_frequencies(f)
S_in = get_positive_frequencies(S_in)
S_out = get_positive_frequencies(S_out)

# Creating some custom subplots
fg, ax = plt.subplots(2,1, figsize=(15,8), sharex=True)

# Plot the absolute values of transformed_output
ax[0].plot(f, np.abs(S_in), 'b-', label='input signal')
ax[0].plot(f, np.abs(S_out), 'r-', label='output1 signal')
ax[0].legend()

ax[1].plot(f, np.abs(S_out/S_in), 'k:', label='transfer function (output/input)')
index=np.argmin(np.abs(np.abs(S_out/S_in)-1/np.sqrt(2)))
plt.axvline(x=index/10,linestyle='--')
print("f0 = {}".format(index/10))
print("inductance = {}".format(200/2/np.pi/index*10))
plt.axhline(y=1/np.sqrt(2),linestyle='--')
ax[1].set_xlabel('Frequency [Hz]')
ax[1].set_ylabel('|T| [V/V]')
ax[1].legend()

ax[0].set_title('FFT of in- and output signals')
ax[0].set_ylabel('|Amplitude| [V/Hz]')

# Fine-tune figure; make subplots close to each other
# and hide x ticks for top plot.
fg.subplots_adjust(hspace=0)
plt.setp(ax[0].get_xticklabels(), visible=False)
plt.show()

f0 = 1328.7
inductance = 0.023956490267463736

mark

1
2
3
# Your solutions to Task 2 here
corner_frequency = 1328.7 # Corner frequency in Hz
inductance = 0.023956490267463736 # Inductance in Henrys

Resonance in RLC series circuits

The sharp peak in transmission for an RLC series circuit associated with a minimum of impedance is referred to as series resonance, with the resonance frequency $f_0$ being that at which the transfer function attains a maximum value (of value 1). With

$$|T(\omega) | = \frac{R}{\sqrt{R^2+(\omega L - 1/\omega C)^2}}$$

this occurs when $\omega L = \frac{1}{\omega C}$ (the inductive reactance equalling the capacitive reactance) and hence

$$\omega_0^2 = \frac{1}{LC} \implies f_0 = \frac{1}{2\pi \sqrt{LC}}$$

For all intents and purposes, a signal is passed by an RLC circuit if it has a frequency contained between the two cut-off frequencies, $f_l$ and $f_h$, at which the transfer function has a magnitude of $\frac{1}{\sqrt{2}}$. The bandwidth of the circuit (also known as the resonance width) is then defined as the difference between these two frequencies, i.e $\ f_h - f_l$.

Task 3

Marks available: 8

You are to work with t3_box in this task (defined in the first cell of the notebook). This contains an RLC band pass filter circuit, with values of resistance, inductance and capacitance different from those of the RLC circuit in Task 1.

Firstly, determine the transfer function for this circuit and plot its magnitude, $|T(f)|$, against frequency, $f$ (in Hz), on the interval $f \in [10,2000]$ Hz on the figure object provided (student_figure). You should ensure that your input signal contains frequency components ranging from 10 to 2000 Hz (remember that you need at least 4000 samples for the frequency domain data to achieve this).

You should use a linear scale on the y-axis (do not use dB), and a logarithmic scale on the frequency axis. This is achieved through the use of the plt.semilogx() function. Your plot should also have a title and axis labels.

Using your plotted data or otherwise, estimate the resonant frequency and the bandwidth of the circuit in Hertz and assign these to the variables resonant_frequency and bandwidth. As in Task 2, these values should be accurate to 10 Hz.

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
# Please do not create any other figure object 
student_figure = plt.figure(figsize=(15,7))

# You may perform additional working, but ensure that the plot
# you wish to be marked is on student_figure defined above

# generate time array
t = get_time_array(N=20000, duration = 10.0)

# generate array containing input signal (voltage)
s_in = swept_sine(10,2000,t)

# Pass input signal through the example_box
# and create an output signal (voltage):
s_out = t3_box.process(t, s_in)

# Perform the Fourier Tranform on both signals
S_in = np.fft.fft(s_in)
S_out = np.fft.fft(s_out)

f = get_frequency_array(t)

f = get_positive_frequencies(f)
S_in = get_positive_frequencies(S_in)
S_out = get_positive_frequencies(S_out)

plt.subplot(211)
plt.plot(f, np.abs(S_in), 'b-', label='input signal')
plt.plot(f, np.abs(S_out), 'r-', label='output1 signal')
plt.title('FFT of in- and output signals')
plt.xlabel(r'Frequency [Hz]')
plt.ylabel(r'|Amplitude| [V/Hz]')

plt.subplot(212)
plt.plot(f, np.abs(S_out/S_in), 'k:', label='transfer function (output/input)')
a=np.abs(S_out/S_in)
fl=np.argmin(np.abs(a[:3000]-1/np.sqrt(2)))/10
fh=np.argmin(np.abs(a-1/np.sqrt(2)))/10
f0=np.argmax(a)
peak=a[f0]
f0/=10
plt.axvline(x=fl,linestyle='--')
plt.axvline(x=fh,linestyle='--')
plt.axvline(x=f0,linestyle='--')
print("fl = {}".format(fl))
print("fh = {}".format(fh))
print("f0 = {}".format(f0))
plt.axhline(y=1/np.sqrt(2),linestyle='--')
plt.axhline(y=peak,linestyle='--')
plt.title('Input Signal')
plt.xlabel(r'Frequency [Hz]')
plt.ylabel(r'|T| [V/V]')

# Automatically adjust subplots to fit figure area
plt.tight_layout()
# Display plots
plt.show()

Data points = 20000
Total time = 9.9995 s
Timestep = 0.0005 s
Sampling Frequency = 2000.0 Hz
fl = 193.3
fh = 614.6
f0 = 344.7

mark

1
2
resonant_frequency = f0 # Resonant frequency in Hz
bandwidth = fh-fl # Bandwidth in Hz

Task 4

Marks available: 5

For the final task you are to create your own black box function called my_box_process(t, s_in, R, L) that mimics the behaviour of an RL low-pass filter circuit of resistance R (Ohms) and inductance L (Henrys).

t and s_in will be two numpy arrays representing an input voltage signal with physical amplitude s_in[i] at time t[i], and the function should return a numpy array containing the real output signal in the time domain.

You will need the transfer function for the RL low-pass filter circuit:
$$T(f) = \frac{1}{1+(f/f_0)j} $$ where $f_0 = R/(2\pi\,L)$.

Hints:

  • Use .fft to perform an FFT on the input signal before applying the transfer function in the frequency domain. You will then need to use the .ifft function to represent this modulated signal in the time domain (ifft in the numpy documentation).
  • You may wish to define the transfer function separately and call this within my_box_process.
  • In python, 1j is used to represent the imaginary unit.
  • Make sure your function returns a physical time domain signal (i.e. real values that can change sign). The ifft function will return an array of complex numbers with the real part of these numbers containing the signal we are interested in.
  • BTW because the input and output of your Box must be numpy arrays, you do not need to use any for loop in your code.
1
2
3
4
5
6
7
8
9
10
11
12
13
# Function definition for Task 4

def my_box_process(t, s_in, R, L):
# Perform the Fourier Tranform on input signals\n",
S_in = np.fft.fft(s_in

f0 = R/2/np.pi/L

T_f = 1/(1+(S_in/f0)*1j

s = np.fft.ifft(T_f)

return s.real
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
# Feel free to insert additional cells for working

# Please do not create any other figure object
student_figure = plt.figure(figsize=(15,7))

# You may perform additional working, but ensure that the plot
# you wish to be marked is on student_figure defined above

# generate time array
t = get_time_array(N=20000, duration = 10.0)

# generate array containing input signal (voltage)
s_in = swept_sine(10,2000,t)

# Pass input signal through the example_box
# and create an output signal (voltage):
s_out = my_box_process(t, s_in, 200, 0.02)

# Perform the Fourier Tranform on both signals
S_in = np.fft.fft(s_in)
S_out = np.fft.fft(s_out)

f = get_frequency_array(t)

f = get_positive_frequencies(f)
S_in = get_positive_frequencies(S_in)
S_out = get_positive_frequencies(S_out)

plt.subplot(211)
plt.plot(f, np.abs(S_in), 'b-', label='input signal')
plt.plot(f, np.abs(S_out), 'r-', label='output1 signal')
plt.title('FFT of in- and output signals')
plt.xlabel(r'Frequency [Hz]')
plt.ylabel(r'|Amplitude| [V/Hz]')

plt.subplot(212)
plt.plot(f, np.abs(S_out/S_in), 'k:', label='transfer function (output/input)')
a=np.abs(S_out/S_in)
fl=np.argmin(np.abs(a[:3000]-1/np.sqrt(2)))/10
fh=np.argmin(np.abs(a-1/np.sqrt(2)))/10
f0=np.argmax(a)
peak=a[f0]
f0/=10
plt.axvline(x=fl,linestyle='--')
plt.axvline(x=fh,linestyle='--')
plt.axvline(x=f0,linestyle='--')
print("fl = {}".format(fl))
print("fh = {}".format(fh))
print("f0 = {}".format(f0))
plt.axhline(y=1/np.sqrt(2),linestyle='--')
plt.axhline(y=peak,linestyle='--')
plt.title('Input Signal')
plt.xlabel(r'Frequency [Hz]')
plt.ylabel(r'|T| [V/V]')

# Automatically adjust subplots to fit figure area
plt.tight_layout()
# Display plots
plt.show()

N.B. Please execute the following statement before submitting your work. It will check that you have used correct naming for the variables and functions specified in the above tasks. It will not tell you whether you have correctly defined and implemented these!

You may execute this statement as many times as required.

1
Boxes.check()

my_box_process function is correctly named.
RL_circuit variable is correctly named.
RC_circuit variable is correctly named.
RLC_circuit variable is correctly named.
corner_frequency variable is correctly named.
inductance variable is correctly named.
resonant_frequency variable is correctly named.
bandwidth variable is correctly named.
RL_circuit, RC_circuit and RLC_circuit are assigned in a valid way.




Interacting with Black Boxes in Assignment 2

Interacting with Black Boxes in Assignment 2

This notebook demonstrates how to interact with the ‘Boxes‘ objects that you will work with in assignment 2. It also includes examples of using the .fft and .fftfreq functions. See also the document Yr2_computing_System_Identification.pdf on Canvas.

Start by executing the code cell below which loads the example Boxes object example_box that represents some unknown electronic circuit.

1
2
3
4
5
6
7
8
9
10
11
# Doing all the imports in one go
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Setting global plotting settings
plt.rcParams['font.size'] = 20

# Creating an example `Box' object to play with
# `Boxes' is part of the assignment module provided
from module_engine.assignment import Boxes
example_box = Boxes.get_example_box()

Using temporary ID. Switch to your own ID as soon as possible.

Generating an input signal

The box object represents a ‘black box’, i.e. an unknown system. The system accepts an input signal and provides an output signal in return. By generating custom input signals for a specific box and analysing the output you can investigate the unknown system.

For the purposes of this example we generate a signal representing an input voltage over time. The signal is generated as a sine wave whose frequency is slowly changed from f1 to f2. In assignment 2, you should consider creating similar signals to probe the behaviour of each box.

1
2
3
4
5
6
7
8
9
10
def swept_sine(f1, f2, t):
"""
create a sine wave with a frequency sweeping between f1 and f2
f1: start frequency [Hz]
f2: stop frequency [Hz]
t: np.array of linearly spaced points in time [s]
"""
f = np.linspace(f1,f2,t.size)
output = np.sin(2*np.pi*f*t)
return output

In order to generate the input signal we first make an array of times (t) at which the signal should be created.

We want to sample the signal on the interval $[0,10)$ seconds with $N = 1000$ sample points, resulting in a spacing between sample points of
timestep $= \frac{10}{1000} = 0.01$ s.

If you think about this in the frequency domain, this choice corresponds to a maximum sampling frequency of $\frac{1}{0.01} = 100$ Hz (i.e. 100 samples per second).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Sampling a signal correctly is a bit tricky, you need
# to make sure you know the total time, timestep and number
# of points accuractely. If in doubt, compute these from the
# time array you created and print the result.

def get_time_array(N=1000, duration=10, endpoint=False):
# N: Number of sample points
# duration: end point of the open interval

t = np.linspace(0, duration, num=N, endpoint=endpoint)
total_time = t[-1]-t[0] # different from duration!!
timestep = total_time/(N-1)
print("Data points = {}".format(N))
print("Total time = {} s".format(total_time))
print("Timestep = {} s".format(timestep))
print("Sampling Frequency = {} Hz".format((1/timestep)))

return t

# generate time array
t = get_time_array(N=1000, duration = 10.0)

# generate array containing input signal (voltage)
s_in = swept_sine(1,10,t)

Data points = 1000
Total time = 9.99 s
Timestep = 0.01 s
Sampling Frequency = 100.0 Hz

Note the endpoint=False keyword parameter. Without it, np.linspace includes the endpoint which would give the unexpected spacing of $1000/999 = 0.01001…$. (The spacing of any numpy array can also be checked using the numpy.diff function: numpy.diff(t) returns an array containing the difference between successive elements of t.)

We now have an input signal. Let’s have a look at it.

1
2
3
4
5
6
7
8
9
fig = plt.figure(figsize=(15,3))

# Plot input voltage as a function of time
plt.plot(t, s_in, 'b-')
plt.title('Input Signal')
plt.xlabel(r'$t$')
plt.ylabel(r'$V_{in}$')
plt.axis([0, 5, np.amin(s_in), np.amax(s_in)])
plt.show()

mark

Probing the box

In the cell, the call example_box.process(t, s_in) passes our voltage signal through the first box and returns the corresponding output signal, stored in the variable s_out. It is essential that the sample points described in t are equally spaced.

Both arrays are then plotted on the interval $[0, 5]$ in the time domain.

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
# Pass input signal through the example_box
# and create an output signal (voltage):
s_out = example_box.process(t, s_in)

# And plot the result, comparing input to output
fig = plt.figure(figsize=(15,8))
# Plot input voltage as a function of time
plt.subplot(211)
plt.plot(t, s_in, 'b-')
plt.title('Input Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{in}$', fontsize=14)
plt.axis([0, 5, np.amin(s_in), np.amax(s_in)])

# Plot output voltage as a function of time
plt.subplot(212)
plt.plot(t, s_out, 'r-')
plt.title('Output Signal', fontsize=14)
plt.xlabel(r'$t$', fontsize=14)
plt.ylabel(r'$V_{out}$', fontsize=14)
plt.axis([0, 5, np.amin(s_out), np.amax(s_out)])

# Automatically adjust subplots to fit figure area
plt.tight_layout()
# Display plots
plt.show()

mark

The signal has certainly been modified by example_box, but it is not yet entirely clear how. This type of investigation is often simpler and more intuitive when describing the signals as a function of frequency, i.e. by using the frequency domain.

To do so, we perform an Fourier Transform (using a DFT algorithm) on both the input and output signals and plot the resulting spectra.

As described in Yr2_Computing_System_Identification.pdf, the function .fftfreq is used to determine the relevant sample frequencies, and we choose to plot the positive frequency terms only.

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
# Perform the Fourier Tranform on both signals
S_in = np.fft.fft(s_in)
S_out = np.fft.fft(s_out)

# Next we generate sample points in the frequency domain,
# for this we need to know the duration and time step
# of the time array.
def get_frequency_array(t):
N = len(t)
timestep = np.diff(t)[0]
f = np.fft.fftfreq(N, d=timestep)
return f

f = get_frequency_array(t)

# We only want to plot the positive frequencies from the FFT
# this means we can ignore half of each array.
# Check the numpy.fftfreq documentation to see that
# frequencies are arranged as follows:
# f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even
# so we are interested in the elements f[1:n/2].
def get_positive_frequencies(A):
N=len(A)
if len(A)//2==0:
N2=int(N/2)
else:
N2=int((N+1)/2)
return A[1:N2]

f = get_positive_frequencies(f)
S_in = get_positive_frequencies(S_in)
S_out = get_positive_frequencies(S_out)

# Creating some custom subplots
fg, ax = plt.subplots(2,1, figsize=(15,8), sharex=True)

# Plot the absolute values of transformed_output
ax[0].plot(f, np.abs(S_in), 'b-', label='input signal')
ax[0].plot(f, np.abs(S_out), 'r-', label='output signal')
ax[0].legend()

ax[1].plot(f, np.abs(S_out/S_in), 'k:', label='transfer function (output/input)')
ax[1].set_xlabel('Frequency [Hz]')
ax[1].set_ylabel('|T| [V/V]')
ax[1].legend()

ax[0].set_title('FFT of in- and output signals')
ax[0].set_ylabel('|Amplitude| [V/Hz]')

# Fine-tune figure; make subplots close to each other
# and hide x ticks for top plot.
fg.subplots_adjust(hspace=0)
plt.setp(ax[0].get_xticklabels(), visible=False)
plt.show()

mark

By dividing the output by the input we can generate very clean data showing the effect of the box independent of the details of the input signal we used. This ratio is called the transfer function of the system, sometimes also called frequency response.

We can see that the box had an effect similar to that of an ideal band pass filter: Our original signal contains frequencies ranging from $1$ to $20$ Hz (see S_in), while the output signal only those from $5$ to $15$ Hz; frequencies outside this range are completely blocked.

You might be curious why our input signal contains significant amplitudes up to frequencies of $20$ Hz even though we specified f2=$10$ Hz. This is because our swept_sine function is just a quick and simple way of creating such a signal. We could instead use a true linear chirp, see the example at the end of this notebook. However, for assignment 2 the exact shape of the signal does not matter.

Additional Notes:

  • In order to complete the first task of assignment 2, you will want to perform the Fourier Transform on the input and output signals as above and calculate the transfer function as their element-wise ratio, .fft(output_signal)/.fft(input_signal) in order to determine the type of circuit contained inside each box.
  • As stated earlier, be sure to investigate a broad range of frequencies with your input signal: 10 - 4000 Hz should be sufficient for Task 1. Additionally, check that your sample interval and number of sample points ($N$) are such that the sample interval in the frequency domains covers all of these frequencies.
  • While the above plots use linear scales, in some cases logarithmic scaling is preferable. For circuits it is also common to plot the transfer function in decibels with a logarithmic horizontal axis. In general: if you are not able to discern anything meaningful from a plot, try using a different axis scaling. For assignment 2, read the instructions carefully and make sure you use the required axis scaling and formatting.

For the curious: how to make a linear chirp

For the tasks in assignment 2 the exact shape of the input signal is not important as long as it contains non-zero amplitudes in the frequency range of interest. The function swept_sine above is a simple way of creating such an input signal. You will have noted above that the output signal of that function shows non-zero amplitudes up to 2 times the maximum frequency (as shown by the Fourier Transform). If you want to create a true linear chirp signal whose frequency sweeps linearly from $f_1$ to $f_2$ you can use the following function. More information on chirp signals can be found online but the short explanation is that the argument of the sine function is the phase $\phi$ and the frequency of a sine is defined as $1/2\pi\,df/d\phi$. Therefore, to create a frequency that increases linear in time, the phase $\phi$ must be a quadratic function of time.

1
2
3
4
5
6
7
8
9
10
def linear_chirp(f1, f2, t):
"""
create a linear chirp with a frequency sweeping between f1 and f2
f1: start frequency [Hz]
f2: stop frequency [Hz]
t: np.array of linearly spaced points in time [s]
"""
k=(f2-f1)/(t[-1]-t[0])
output = np.sin(2*np.pi*(f1*t + k/2.0*t*t))
return output

Let’s try it out:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
t = get_time_array(N = 500, duration = 10)
s_in = linear_chirp(1,10,t)

f = get_frequency_array(t)
S_in = np.fft.fft(s_in)

f = get_positive_frequencies(f)
S_in = get_positive_frequencies(S_in)

fig = plt.figure(figsize=(15,3))
plt.plot(f, np.abs(S_in), 'b-')
plt.xlabel('Frequency [Hz]')
plt.ylabel('|S_in|')
plt.show()

Data points = 500
Total time = 9.98 s
Timestep = 0.02 s
Sampling Frequency = 50.0 Hz

mark

×

纯属好玩

扫码支持
扫码打赏,你说多少就多少

打开支付宝扫一扫,即可进行扫码打赏哦

文章目录
  1. 1. ASSIGNMENT 2:System Identification
    1. 1.1. Task 1
      1. 1.1.1. Marks available: 4
    2. 1.2. Task 2
      1. 1.2.1. Marks available: 3
    3. 1.3. Resonance in RLC series circuits
    4. 1.4. Task 3
      1. 1.4.1. Marks available: 8
    5. 1.5. Task 4
      1. 1.5.1. Marks available: 5
  2. 2. Interacting with Black Boxes in Assignment 2
    1. 2.1. Interacting with Black Boxes in Assignment 2
      1. 2.1.1. Generating an input signal
      2. 2.1.2. Probing the box
    2. 2.2. For the curious: how to make a linear chirp