import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import chirp, max_len_seq, freqz, fftconvolve, resample
import sounddevice as sd
import tools
def crest_factor(x):
"""Peak-to-RMS value (crest factor) of the signal x
Parameter
---------
x : array_like
signal
"""
return np.max(np.abs(x)) / np.sqrt(np.mean(x**2))
def circular_convolve(x, y, outlen):
"""Circular convolution of x and y
Parameters
----------
x : array_like
Real-valued signal
y : array_like
Real-valued signal
outlen : int
Length of the output
"""
return np.fft.irfft(np.fft.rfft(x, n=outlen) * np.fft.rfft(y, n=outlen), n=outlen)
def plot_time_domain(x, fs=44100, ms=False):
time = np.arange(len(x)) / fs
timeunit = 's'
if ms:
time *= 1000
timeunit = 'ms'
fig = plt.figure()
plt.plot(time, x)
plt.xlabel('Time / {}'.format(timeunit))
return
def plot_freq_domain(x, fs=44100, khz=False):
Nf = len(x) // 2 + 1
freq = np.arange(Nf) / Nf * fs / 2
frequnit = 'Hz'
if khz:
freq /= 1000
frequnit = 'kHz'
fig = plt.figure()
plt.plot(freq, db(np.fft.rfft(x)))
plt.xscale('log')
plt.xlabel('Frequency / {}'.format(frequnit))
plt.ylabel('Magnitude / dB')
return
def compare_irs(h1, h2, ms=False):
t1 = np.arange(len(h1)) / fs
t2 = np.arange(len(h2)) / fs
timeunit = 's'
if ms:
t1 *= 1000
t2 *= 1000
timeunit = 'ms'
fig = plt.figure()
plt.plot(t1, h1, t2, h2)
plt.xlabel('Time / {}'.format(timeunit))
return
def compare_tfs(h1, h2, khz=False):
n1 = len(h1) // 2 + 1
n2 = len(h2) // 2 + 1
f1 = np.arange(n1) / n1 * fs / 2
f2 = np.arange(n2) / n2 * fs / 2
frequnit = 'Hz'
if khz:
freq /= 1000
frequnit = 'khz'
fig = plt.figure()
plt.plot(f1, db(np.fft.rfft(h1)), f2, db(np.fft.rfft(h2)))
plt.xscale('log')
plt.xlabel('Frequency / {}'.format(frequnit))
plt.ylabel('Magnitude / dB')
return
def pad_zeros(x, nzeros):
"""Append zeros at the end of the input sequence
"""
return np.pad(x, (0, nzeros), mode='constant', constant_values=0)
fs = 44100
dur = 1
L = int(np.ceil(dur * fs))
time = np.arange(L) / fs
Generate a random signal with normal (Gaussian) amplitude distribution. Use numpy.random.randn
and normalize the amplitude with tools.normalize
.
Let's listen to it.
Plot the signal in the time domain and in the frequency domain.
Is the signal really white?
What is the crest factor of a white noise?
Now feed the white noise to an unkown system tools.blackbox
and save the output signal.
How do you think we can extract the impulse response of the system?
Try to compute the impulse response from the output signal.
Compare it with the actual impulse response which can be obtained by feeding an ideal impulse to tools.blackbox
.
Maximum-length sequences (MLSs) are binary sequences that can be generated very easily with an N-staged shift register and an XOR gate (with up to four inputs) connected with the shift register in such a way that all possible 2N states, minus the case "all 0," are run through. This can be accomplished by hardware with very few simple TTL ICs or by software with less than 20 lines of assembly code.
(Müller 2001)
nbit = int(np.ceil(np.log2(L)))
mls, _ = max_len_seq(nbit) # sequence of 0 and 1
mls = 2*mls - 1 # sequence of -1 and 1
Take a look at the signal in the time domain.
Examine the properties of the MLS
tools.blackbox
In practice, the (digital) signal has to be converted into an analog signal by an audio interface? Here, the process is simulated by oversampling the signal by a factor of 10. Pay attention to the crest factor before and after upsampling.
upsample = 10
mls_up = resample(mls, num=len(mls) * upsample)
time = np.arange(len(mls)) / fs
time_up = np.arange(len(mls_up)) / fs / upsample
plt.figure(figsize=(10, 4))
plt.plot(time_up, mls_up, '-', label='Analog')
plt.plot(time, mls, '-', label='Digital')
plt.legend(loc='best')
plt.xlabel('Time / s')
plt.title('Crest factor {:.1f} -> {:.1f} dB'.format(
tools.db(crest_factor(mls)), tools.db(crest_factor(mls_up))))
plt.figure(figsize=(10, 4))
plt.plot(time_up, mls_up, '-', label='Analog')
plt.plot(time, mls, 'o', label='Ditigal')
plt.xlim(0, 0.0025)
plt.legend(loc='best')
plt.xlabel('Time / s')
plt.title('Crest factor {:.1f} -> {:.1f} dB'.format(
tools.db(crest_factor(mls)), tools.db(crest_factor(mls_up))))
Generate a linear sweep with lin_sweep
.
def lin_sweep(fstart, fstop, duration, fs):
"""Generation of a linear sweep signal.
Parameters
----------
fstart : int
Start frequency in Hz
fstop : int
Stop frequency in Hz
duration : float
Total length of signal in s
fs : int
Sampling frequency in Hz
Returns
-------
array_like
generated signal vector
Note that the stop frequency must not be greater than half the
sampling frequency (Nyquist-Shannon sampling theorem).
"""
if fstop > fs / 2:
raise ValueError("fstop must not be greater than fs/2")
t = np.arange(0, duration, 1 / fs)
excitation = np.sin(
2 * np.pi * ((fstop - fstart) /
(2 * duration) * t ** 2 + fstart * t))
# excitation = excitation - np.mean(excitation) # remove direct component
return excitation
fs = 44100
fstart =
fstop =
duration =
lsweep =
Examine the properties of linear sweeps
pyplot.specgram
with NFFT=512
and Fs=44100
)tools.blackbox
Generate a exponential sweep with exp_sweep
.
def exp_sweep(fstart, fstop, duration, fs):
"""Generation of a exponential sweep signal.
Parameters
----------
fstart : int
Start frequency in Hz
fstop : int
Stop frequency
duration : float
Total length of signal in s
fs : int
Sampling frequency in Hz
Returns
-------
array_like
Generated signal vector
Note that the stop frequency must not be greater than half the
sampling frequency (Nyquist-Shannon sampling theorem).
"""
if fstop > fs / 2:
raise ValueError("fstop must not be greater than fs/2")
t = np.arange(0, duration, 1 / fs)
excitation = np.sin(2 * np.pi * duration *
fstart / np.log(fstop / fstart) *
(np.exp(t / duration * np.log(fstop / fstart)) - 1))
# excitation = excitation - np.mean(excitation) # remove direct component
return excitation
fs = 44100
fstart =
fstop =
duration =
esweep =
Examine the properties of linear sweeps
pyplot.specgram
with NFFT=512
and Fs=44100
)tools.blackbox