MP1: Signals and Systems

This machine problem will explore the different attributes that a system can have: linearity, time-invariance, causality or anti-causality, stability.

We will look at signals generated by an inertial measurement unit (IMU), specifically we will be looking at the accelerometer signals, which tell us the orientation and acceleration of the device, and the gyroscope signals, which tell us the rotational velocity of the device. If somebody's wearing the device, these things tell us:

• Are they standing, walking, sitting, reclining, or lying down?
• If they're walking, how fast?
• If they're sleeping, in what position are they lying?

The systems that extract such information from IMU signals including a complicated mixture of linear vs. nonlinear, time-varying vs. time-invariant, and causal vs. non-causal systems.

Contents

1. Reading the Data
2. Processing the Data
3. Linearity
4. Time Invariance
5. Causality
6. Stability
7. Submitting your answers to the autograder

# Here are some standard packages that we will use almost all the time
import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
import importlib

1. Reading the Data

Scikit-digital-health (SKDH) is a toolkit with IO functions that allow you to read in data from a number of different types of smartwatches and other IMU devices. You can install it using something like the following:

!pip install -q scikit-digital-health

In this MP we'll use data from the Newcastle Polysomnography and Accelerometer Dataset, which is distributed under a Creative Commons attribution license. NPAD contains recordings of 28 participants, including polysomnography (a device with about two dozen electrodes measuring many physiological functions) and accelerometers worn on both the left and right wrist. The accelerometers were GeneActiv devices, so they can be read using an SKDH ReadBin object with default settings:

import skdh.io
reader = skdh.io.ReadBin()

df = reader.predict(file='MECSLEEP42_left_wrist.bin')
print('df is a dict with the following keys:', list(df.keys()))
print('sampling frequency is:',df['fs'])
print('time is an array containing',len(df['time']),'timestamps')
print('accel is an array of shape',df['accel'].shape)
print('the first five samples of accel are:')
print(df['accel'][:5,:])

df is a dict with the following keys: ['accel', 'temperature', 'light', 'fs', 'time']
sampling frequency is: 85.7
time is an array containing 6333900 timestamps
accel is an array of shape (6333900, 3)
the first five samples of accel are:
[[ 0.21858754 -0.01966837 -0.98852907]
 [ 0.19859234 -0.01571095 -0.99243075]
 [ 0.21058946 -0.01571095 -0.99243075]
 [ 0.21058946 -0.03549804 -0.96902068]
 [ 0.21058946 -0.01571095 -1.0002341 ]]

/opt/anaconda3/lib/python3.13/site-packages/skdh/io/geneactiv.py:106: RuntimeWarning: Block fs is not the same as header fs. Setting to block fs.
  n_max, fs, acc, time, light, temp = read_geneactiv(str(file))
/opt/anaconda3/lib/python3.13/site-packages/skdh/io/base.py:111: UserWarning: Timestamps are local but naive, and no time-zone information is available. This may mean that if a DST change occurs during the recording period, the times will be offset by an hour
  warn(

An accelerometer records three signals in parallel: acceleration in the X, Y, and Z directions. The GeneActiv accelerometer normalizes all measurements so that 1.0 is the acceleration of gravity (9.8 m/s^2). To better understand what's happening, let's plot the data.

axs = plt.subplots(2, figsize=(14,6), layout='tight')
print(axs)
axs[1][0].plot(df['time']-df['time'][0], df['accel'])
axs[1][0].legend(['X','Y','Z'])
axs[1][0].set_title('Accelerometers (g)')
axs[1][1].plot(df['time']-df['time'][0], df['temperature'])
axs[1][1].set_title('Temperature (c)')
axs[1][1].set_xlabel('Time (s)')

(<Figure size 1400x600 with 2 Axes>, array([<Axes: >, <Axes: >], dtype=object))

Text(0.5, 0, 'Time (s)')

/opt/anaconda3/lib/python3.13/site-packages/IPython/core/events.py:82: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  func(*args, **kwargs)
/opt/anaconda3/lib/python3.13/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)

Some low temperatures may suggest that the accelerometer was only worn for part of the day. We can check that using SKDH's implementation of the "DETACH" wear-detect algorithm:

df.update(skdh.preprocessing.DETACH().predict(time=df['time'],accel=df['accel'],temperature=df['temperature'],fs=df['fs']))
print('The accelerometer wear periods are:')
print(df['wear']/df['fs'])

The accelerometer wear periods are:
[[    0.          1405.42590432]
 [ 5776.4294049  71950.46674446]]

Let's plot just the part when the watch is worn. Let's also plot only the non-artifact part of the accelerometer signal, i.e., the part with amplitude less than 2g:

ns = df['wear'][1][0]
ne = df['wear'][1][1]

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((df['time'][ns:ne]-df['time'][0])/3600, df['accel'][ns:ne,:])
ax[1].set_ylim([-2,2])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Accelerometer signals')

Text(0.5, 1.0, 'Accelerometer signals')

The timestamps are in seconds relative to the start of the unix era (January 1, 1970), so they can be used to figure out on which day or days this file was recorded, and divide it into days. In sleep studies, it is standard to start every day at noon (12:00), and set it to last for 24 hours:

df.update(skdh.preprocessing.GetDayWindowIndices(bases=12, periods=24).predict(time=df['time']))
df

{'accel': array([[ 0.21858754, -0.01966837, -0.98852907],
        [ 0.19859234, -0.01571095, -0.99243075],
        [ 0.21058946, -0.01571095, -0.99243075],
        ...,
        [ 0.05462689, -0.97736357, -0.15357004],
        [ 0.00663841, -1.02881   , -0.18868513],
        [-0.033352  , -0.96549131, -0.17307842]], shape=(6333900, 3)),
 'temperature': array([26.9, 26.9, 26.9, ..., 22.9, 22.9, 22.9], shape=(6333900,)),
 'light': array([ 0.,  0.,  0., ..., 85., 85., 88.], shape=(6333900,)),
 'fs': 85.7,
 'time': array([1.4259240e+09, 1.4259240e+09, 1.4259240e+09, ..., 1.4259979e+09,
        1.4259979e+09, 1.4259979e+09], shape=(6333900,)),
 'wear': array([[      0,  120445],
        [ 495040, 6166155]]),
 'day_ends': {(np.int64(12),
   np.int64(24)): array([[      0, 5554243],
         [5554243, 6333899]])}}

We see that this file includes recordings from two separate days: The first 5.6M samples are from one day (until noon), the last 0.7M samples are from the next day (after noon on the next day). To find out exactly which days those were, we can run the Sleep().predict(), which will also tell us when the person fell asleep (TSO), how long they slept (in minutes), the number of times they woke up, the number of minutes slept during each sleep bout, and the number of minutes awake during each wake bout:

proc = skdh.sleep.Sleep()
df2 = proc.predict(time=df['time'],
                   accel=df['accel'],
                   temperature=df['temperature'], 
                   fs=df['fs'], 
                   day_ends=df['day_ends'],
                   wear=df['wear'])
df2

{'Day N': [1],
 'Date': ['2015-03-09'],
 'Day Start Timestamp': ['2015-03-09 18:00:00.500000'],
 'Day End Timestamp': ['2015-03-10 12:00:00.000000'],
 'Total Minutes': [np.float64(1080.0)],
 'Wear Minutes': [np.float64(1007.2)],
 'TSO Start Timestamp': [np.float64(1425930835.5)],
 'TSO Start': ['2015-03-09 19:53:55.500000'],
 'TSO Duration': [np.float64(639.25)],
 'total sleep time': [np.int64(618)],
 'percent time asleep': [np.float64(96.562)],
 'number of wake bouts': [np.int64(10)],
 'sleep onset latency': [np.int64(1)],
 'wake after sleep onset': [np.int64(20)],
 'average sleep duration': [np.float64(56.18181818181818)],
 'average wake duration': [np.float64(1.8333333333333333)],
 'sleep wake transition probability': [np.float64(0.01779935275080906)],
 'wake sleep transition probability': [np.float64(0.5454545454545455)],
 'sleep gini index': [np.float64(0.6776699029126213)],
 'wake gini index': [np.float64(0.40495867768595045)],
 'sleep average hazard': [np.float64(0.3322871572871573)],
 'wake average hazard': [np.float64(0.8055555555555555)],
 'sleep power law distribution': [np.float64(1.281471126954279)],
 'wake power law distribution': [np.float64(1.9617966939259757)]}

2. Processing the Data

As suggested above, skdh.sleep.Sleep().predict() estimates the wearer's level of physical activity every minute, thresholds the activity levels, then groups similar regions to identify possible sleep bouts. The algorithms for doing so are interesting; if you want to learn more about them, you can read the code. In this lab we will do all of those same steps, but using the algorithm from a different paper. We will follow the LIDS system (Locomotor Inactivity During Sleep) as reported in the paper "Dynamics and Ultradian Structure of Human Sleep in Real Life" (Eva Winnebeck, Dorothee Fischer, Tanya Leise & Till Roenneberg, 2018).

It is useful to line up the two algorithms, to see what they have in common. Both algorithms include the following steps:

1. Lowpass filter the signal to eliminate high-frequency noise
2. Downsample (SKDH: to 20Hz, LIDS: to 1Hz)
3. Highpass filter (SKDH) or subtract the mean (LIDS)
4. Activity level = RMS magnitude of the acccelerometer vector
5. Compute a local 7-minute weighted local sum (SKDH) or 10-minute unweighted local sum (LIDS)
6. Propose a sleep bout whenever the local summed activity is below 0.5g (SKDH) or below 15% of the full-day average (LIDS)
7. Require that a sleep bout must be at least 30m (LIDS) or that a wake bout after 15m of sleep must be at least 5m (SKDH)
8. LIDS only: LIDS = 100/(activity level + 1)

flowchart LR
    X[(accel)] --> A
    A[LPF] --> B[DS]
    B --> C[HPF]
    C --> D[RMS]
    D --> E[LocSum]
    E --> G[Threshold]
    G --> H[MinDur]
    H --> I[Invert]
    I --> Y[(LIDS)]

This section of the Jupyter notebook will define methods for each of the blocks in the flowchart above. Your homework, discussed in the next four sections of this notebook, will be to create functions testing the linearity, time invariance, causality, and stability of each of those functions.

2.1 LPF

A lowpass filter smooths the input waveform, eliminating high-frequency noise. An LPF can be implemented in many different ways, some of which you'll learn during the remainder of this course. For now, let's use scipy.signal.butter to design the filter, and scipy.signal.sosfilter to implement it.

def section1_lpf(lpf_input, cutoff_frequency=1/df['fs'], axis=0):
    sos = scipy.signal.butter(3, cutoff_frequency, btype='low', output='sos')
    lpf_output = scipy.signal.sosfilt(sos, lpf_input, axis)
    return lpf_output

For successful downsampling (as you will learn next week!), we need to lowpass filter with a cutoff frequency equal to one half-cycle per $W$ samples, where

$$W=\frac{F_{s,output}}{F_{s,input}}$$

is the ratio of the output sampling rate divided by the input sampling rate. In our case, $F_{s,output}=1$, and $F_{s,input}=$df['fs'].

lpf_output = section1_lpf(df['accel'])

ns = df['wear'][1][0]
ne = df['wear'][1][1]

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((df['time'][ns:ne]-df['time'][0])/3600, lpf_output[ns:ne,:])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Lowpass-filtered accelerometer signals')

Text(0.5, 1.0, 'Lowpass-filtered accelerometer signals')

2.2 Downsample

A lowpass filter followed by downsampling is called decimation. We could have combined steps 1 and 2 by using scipy.signal.decimate, but we will keep them separate in this notebook so that we can separately analyze their system properties.
Downsampling just means that we throw away $D-1$ out of every $D$ samples, where

$$D=\frac{1}{W}=\frac{F_{s,input}}{F_{s,output}}$$

If $D$ is not an integer, we can get better results by resampling using scipy.signal.resample or scipy.signal.resample_poly. In this notebook we will use an approximate method in which we keep $x[n]$ only if $n//D\ne (n-1)//D$, where $//$ is integer division. This method guarantees that the output signal has the correct sampling rate, but it means that the time alignment of the samples is a little bit different from one sample to the next.

downsampled_fs = 1

def section2_downsample(downsample_input, D=df['fs']/downsampled_fs):
    downsample_output = [ downsample_input[0] ]
    for n in range(1,len(downsample_input)):
        if n//D != (n-1)//D:
            downsample_output.append(downsample_input[n])
    return np.array(downsample_output)

downsample_output = section2_downsample(lpf_output)
downsample_time = section2_downsample(df['time'])

print('The shape of the signal before downsampling is:',lpf_output.shape)
print('The shape of the signal after downsampling is:',downsample_output.shape)

The shape of the signal before downsampling is: (6333900, 3)
The shape of the signal after downsampling is: (73908, 3)

ns = int(df['wear'][1][0]*downsampled_fs/df['fs'])
ne = int(df['wear'][1][1]*downsampled_fs/df['fs'])

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((downsample_time[ns:ne]-downsample_time[0])/3600, downsample_output[ns:ne,:])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Accelerometer signals downsampled to 1 sample/second')

Text(0.5, 1.0, 'Accelerometer signals downsampled to 1 sample/second')

2.3 HPF or subtract the local average

Next, we want to emphasize only the local changes in the accelerometer readings. SKDH does this by highpass filtering with a cutoff frequency of 0.25Hz. LIDS gets roughly the same effect by calculating the local average of the signal, in ten-minute chunks, then subtracting the local average from each sample. Let's do the latter.

def section3_hpf(hpf_input, nsamps=600*downsampled_fs, axis=0):
    '''Subtract the local average over a period of nsamps samples'''
    hpf_output = np.zeros(hpf_input.shape)
    for t in range(len(hpf_input)):
        locave = np.average(hpf_input[max(t-nsamps//2,0):min(t-nsamps//2+nsamps,len(hpf_input))], axis=axis)
        hpf_output[t] = hpf_input[t] - locave
    return hpf_output

Ten minutes = 600 seconds times the downsampled sampling rate:

hpf_output = section3_hpf(downsample_output)

print('hpf_output shape is',hpf_output.shape)

hpf_output shape is (73908, 3)

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((downsample_time[ns:ne]-downsample_time[0])/3600, hpf_output[ns:ne,:])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Accelerometer signals minus their 10-minute local averages')

Text(0.5, 1.0, 'Accelerometer signals minus their 10-minute local averages')

2.4 RMS magnitude

"Actigraphy" is the study of human rest and activity cycles. Its primary object of study is the "activity level," which is the amount by which a person moves in any given second, plotted as a function of time.

An accelerometer measures acceleration in three dimensions. The most important acceleration acting on the human body at any given time is gravity, therefore what an accelerometer really measures is orientation relative to the direction of gravity. By highpass-filtering or subtracting the mean, one obtains a signal showing only the change in orientation from one second to the next. "Activity level" is the RMS magnitude of that change. If $a_x$, $a_y$, and $a_z$ are the X, Y, and Z components of acceleration, with average values of $\bar{a}_x,\bar{a}_y,\bar{a}_z]$, then the activity level is

$$a = |\mathbf{a}-\bar{\mathbf{a}}| = \sqrt{(a_x-\bar{a}_x)^2+(a_y-\bar{a}_y)^2+(a_z-\bar{a}_z)^2}$$

def section4_rms(rms_input, axis=1):
    return np.sqrt(np.sum(np.square(rms_input), axis=axis))

rms_output = section4_rms(hpf_output)

print('rms_output shape is',rms_output.shape)

rms_output shape is (73908,)

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((downsample_time[ns:ne]-downsample_time[0])/3600, rms_output[ns:ne])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Activity level per second')

Text(0.5, 1.0, 'Activity level per second')

2.5 Local Sum

Different aspects of actigraphy look for different patterns in the activity signal:

• Gait analysis finds periodic fluctuations, with a period similar to footsteps, i.e., about one second
• Physical activity analysis looks for regions of low, medium, and high activity
• Sleep analysis looks for prolonged periods of very low activity, separated by brief wake bouts

To find prolonged periods of very low activity, the first thing is to compute a local sum over time windows of about 7-10 minutes (420-600 seconds)

def section5_locsum(locsum_input, nsamps=600*downsampled_fs):
    '''Compute unweighted local sum over nsamps samples'''
    locsum_output = np.zeros(locsum_input.shape)
    for n in range(len(locsum_input)):
        locsum_output[n] = np.sum(locsum_input[max(0,n-nsamps//2):min(len(locsum_output),n-nsamps//2+nsamps)])
    return locsum_output

locsum_output = section5_locsum(rms_output)

print('locsum_output shape is',locsum_output.shape)

locsum_output shape is (73908,)

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((downsample_time[ns:ne]-downsample_time[0])/3600, locsum_output[ns:ne])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('10-minute local summed activity level')

Text(0.5, 1.0, '10-minute local summed activity level')

2.6 Threshold

If a person is sufficiently inactive, for a sufficiently long period of time, we propose that they might be asleep. LIDS proposes that a person might be asleep if their activity level is below 15 percent of the daily average; the next step will set a minimum duration.

def section6_threshold(threshold_input):
    '''Output a 1 if the input is below 15 percent of average, otherwise 0'''
    threshold = 0.15 * np.average(threshold_input)
    threshold_output = np.array([ 1 if threshold_input[n] < threshold else 0 for n in range(len(threshold_input)) ])
    return threshold_output

threshold_output = section6_threshold(locsum_output)

print('threshold_output shape is',threshold_output.shape)

threshold_output shape is (73908,)

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((downsample_time[ns:ne]-downsample_time[0])/3600, threshold_output[ns:ne])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Candidate sleep bouts')

Text(0.5, 1.0, 'Candidate sleep bouts')

2.7 Minimum Duration

Being still for a few minutes does not imply sleep. Maybe we can say that you are asleep if you are still for at least 30 minutes.

def section7_mindur(mindur_input, nsamps=30*60*downsampled_fs):
    '''
    Find the onset and offset of every candidate sleep bout.  
    If offset-onset>nsamps, keep the sleep bout, otherwise delete it.
    '''
    mindur_output = np.zeros(len(mindur_input))
    onsets = np.where(np.diff(mindur_input)>0)[0]
    offsets = np.where(np.diff(mindur_input)<0)[0]
    for k in range(min(len(onsets),len(offsets))):
        if offsets[k]-onsets[k] > nsamps:
            mindur_output[onsets[k]+1:offsets[k]+1] = 1
    return mindur_output

mindur_output = section7_mindur(threshold_output)

print('mindur_output shape is',mindur_output.shape)

mindur_output shape is (73908,)

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((downsample_time[ns:ne]-downsample_time[0])/3600, mindur_output[ns:ne])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Confirmed sleep bouts')

Text(0.5, 1.0, 'Confirmed sleep bouts')

2.8 Regularized percentage inverse

Winnebeck et al. (2018) proposed that it's possible to differentiate REM sleep (rapid eye movement sleep, i.e., "active" dreaming sleep) versus NREM sleep (not REM) by computing the "locomotor inactivity during sleep" (LIDS), which is the regularized percentage inverse of activity during sleep, and then using a Fourier transform to identify periodic variation in this measure. We will learn about Fourier transforms later this semester, but we can learn about regularized percentage inverses right now. The word "regularized" means that you add 1 to the denominator so there are no divide-by-zero errors; the word "percentage" means multiply by 100; the word "inverse" just means one over the activity, so

$$\text{LIDS}[t]=\frac{100}{1+a[t]}$$

where $a[t]$ is the activity (locsum_output).

def section8_lids(activity):
    return 100 / (1+activity)

lids_output = section8_lids(locsum_output)

print('lids_output shape is',lids_output.shape)

lids_output shape is (73908,)

ax = plt.subplots(1, figsize=(14,4), layout='tight')
ax[1].plot((downsample_time[ns:ne]-downsample_time[0])/3600, lids_output[ns:ne])
ax[1].set_xlabel('Time (hours)')
ax[1].set_ylabel('Accel (fractions of g)')
ax[1].set_title('Locomotor inactivity during sleep')

Text(0.5, 1.0, 'Locomotor inactivity during sleep')

Finding REM vs NREM sleep for this person's record is difficult, because they woke up so often during the night. There are little dips in the LIDS measure that might indicate REM sleep; more careful analysis using a Fourier transform would be needed to confirm this hypothesis... But that's not the purpose of this MP. The purpose of this MP is to figure out which of our system components are linear, which are time-invariant, which are causal, and which are stable.

3. Linearity

A system $H$ is linear if and only if, for any two signals $x_1[n]$ and $x_2[n]$ and for any scalar $a$,

$$x_1[n]\rightarrow \boxed{H}\rightarrow y_1[n]$$ $$x_2[n]\rightarrow \boxed{H}\rightarrow y_2[n]$$ $$x_1[n]+ax_2[n]\rightarrow \boxed{H}\rightarrow y_1[n]+ay_2[n]$$

• You can't use software to demonstrate that $H$ is linear, because you can't test it for every possible input
• You can use software to demonstrate that $H$ is not linear: if it fails the test for any signal you can think of, then it's not linear

import submitted
importlib.reload(submitted)
help(submitted.test_linearity)

Help on function test_linearity in module submitted:

test_linearity(system, x1, x2, a, tol=1e-05)
    Run a test, using provided signals, to determine whether the system might be linear.

    @param:
    system - a callable object, e.g., a method
    x1 - an ndarray
    x2 - an ndarray of the same size as x1
    a - a scalar
    tol - consider two signals equal if the average abs difference is less than tol

    @return:
    True if y3=y1+a*y2, else False, where y3=system(x1+a*x2).
    Here "True" means the system might be linear, "False" means it definitely isn't.

Open submitted.py in an editing window, and enter code to implement the test described above. When you are done, you should get results like this:

importlib.reload(submitted)

# Define x1 and x2, kind of arbitrarily, to be the first and second thousand 35000 samples of downsample_output
x1 = downsample_output[:35000,:]  
x2 = downsample_output[35000:70000,:]

print(submitted.test_linearity(section1_lpf, x1, x2, 0.5), ': section1_lpf might be linear.')
print(submitted.test_linearity(section2_downsample, x1, x2, 0.5), ': section2_downsample might be linear')
print(submitted.test_linearity(section3_hpf, x1, x2, 0.5), ': section3_hpf might be linear')
print(submitted.test_linearity(section4_rms, x1, x2, 0.5), ': section4_rms might be linear')
print(submitted.test_linearity(section5_locsum, x1[:,0], x2[:,0], 0.5), ': section5_locsum might be linear')
print(submitted.test_linearity(section6_threshold, x1[:,0], x2[:,0], 0.5), ': section6_threshold might be linear')
print(submitted.test_linearity(section7_mindur, x1[:,0], x2[:,0], 0.5), ': section7_mindur might be linear')
print(submitted.test_linearity(section8_lids, x1[:,0], x2[:,0], 0.5), ': section8_lids might be linear')

True : section1_lpf might be linear.
True : section2_downsample might be linear
True : section3_hpf might be linear
False : section4_rms might be linear
True : section5_locsum might be linear
False : section6_threshold might be linear
True : section7_mindur might be linear
False : section8_lids might be linear

The results are:

• Nonlinear systems: RMS, Threshold, LIDS
• Possibly linear systems: LPF, Downsample, HPF, Local Sum, Minimum Duration

Proving that LPF is linear and time-invariant requires some knowledge of filters, which will be covered later in the course. Proving linearity or nonlinearity of all other systems, however, is just algebra: You should be able to do it already right now!

Can you use pencil and paper to prove that RMS, Threshold, and LIDS are nonlinear? Can you prove that Downsample, HPF, and Local Sum are linear? Is MinDur really linear, or are there signals for which it would not pass the linearity test?

4. Time Invariance

A system $H$ is time invariant if and only if, for any signal $x[n]$ and integer $d$,

$$x[n]\rightarrow \boxed{H}\rightarrow y[n]$$ $$x[n-d]\rightarrow \boxed{H}\rightarrow y[n-d]$$

Usually, you don't use software to demonstrate that the equations above are true for any $x[n]$ and $d$. Instead, you can test linearity for one particular $x[n]$ and $d$. If the test fails, you know the system is not time invariant. If the test succeeds, you don't know that the system is time invariant, but you at least know if behaves time-invariant for the signal that you provided.

We will pass the "system" in the form of a python method. Load submitted.py, and read the docstring for test_timeinvariance:

importlib.reload(submitted)
help(submitted.test_timeinvariance)

Help on function test_timeinvariance in module submitted:

test_timeinvariance(system, x, d, tol=1e-05)
    Run a test, using the provided signal, to determine whether the system might be time-invariant:
    Create a signal x2 that has the same shape as x, but x2[n]=x[n-d] for d<=n<len(x).
    Create y2=system(x2) and y=system(y), then test to see if y2[n]=y[n-d] for d<=n<len(y).
    If not, the system is not time-invariant; if so, it might be.

    @param:
    system - a callable object, e.g., a method
    x - an ndarray
    d - a positive integer
    tol - consider two signals equal if the average abs difference is less than tol

    @return:
    True if the system might be time-invariant, otherwise False.

The test described by this docstring has a bug: It may behave strangely if the signal has important components near the beginning or near the end. In order to avoid such problems, let's use a test signal composed of 180,000 rows of zeros, followed by the signal of interest, followed by another 180,000 rows of zeros.

Open submitted.py in an editing window, and enter code to implement the test described above. When you are done, you should get results like this:

importlib.reload(submitted)

# Define x to equal downsample_input, preceded and followed by 180k rows of zeros
x = np.zeros((360000+len(downsample_output), 3))
x[180000:180000+len(downsample_output),:] = downsample_output

print(submitted.test_timeinvariance(section1_lpf, x, 50), ': section1_lpf might be time invariant.')
print(submitted.test_timeinvariance(section2_downsample, x, 50), ': section2_downsample might be time invariant')
print(submitted.test_timeinvariance(section3_hpf, x, 50), ': section3_hpf might be time invariant')
print(submitted.test_timeinvariance(section4_rms, x, 50), ': section4_rms might be time invariant')
print(submitted.test_timeinvariance(section5_locsum, x[:,0], 50), ': section5_locsum might be time invariant')
print(submitted.test_timeinvariance(section6_threshold, x[:,0], 50), ': section6_threshold might be time invariant')
print(submitted.test_timeinvariance(section7_mindur, x[:,0], 50), ': section7_mindur might be time invariant')
print(submitted.test_timeinvariance(section8_lids, x[:,0], 50), ': section8_lids might be time invariant')

True : section1_lpf might be time invariant.
False : section2_downsample might be time invariant
True : section3_hpf might be time invariant
True : section4_rms might be time invariant
True : section5_locsum might be time invariant
True : section6_threshold might be time invariant
True : section7_mindur might be time invariant
True : section8_lids might be time invariant

The results are:

• Downsampling is not time-invariant
• All of the other systems might be time-invariant!

Can you use pencil and paper to prove that downsampling is time-varying? Can you use pencil and paper to prove that HPF, RMS, Local Sum, Threshold, Min Duration, and LIDS are time-invariant? LPF is also time-invariant, but proving that it's time-invariant (for every possible input) requires knowledge of filtering that you will learn later in the course.

5. Causality

A system is causal if and only if the output, $y[n]$, depends on only current and past values of the input.

If the input and output have the same sampling rate, the definition above can be written as: $y[n]$ depends on values of $x[m]$ only for $m\le n$. All of our systems have this property except Downsampling. Downsampling is also causal, but we won't prove it here.

importlib.reload(submitted)
help(submitted.test_causality)

Help on function test_causality in module submitted:

test_causality(system, x, n, tol=1e-05)
    Run a test, using the provided signal, to determine whether the system might be causal:
    Check to see whether y[n]=y2[n], where
     - y[n] is generated using the entire x, including past and future parts
     - y2[n] is generated using only x[m] for m<=n.
    If y[n]=y2[n], then the system might be causal; if not, it's non-causal.

    @param:
    system - a callable object, e.g., a method
    x - an ndarray
    n - a positive integer
    tol - consider y[n] and y2[n] equal if the np.linalg.norm(y[n]-y2[n]) < tol

    @return:
    True if the system might be causal, False if it's konwn to be non-causal.

importlib.reload(submitted)

x = downsample_output

print(submitted.test_causality(section1_lpf, x, 500), ': section1_lpf might be causal.')
print(submitted.test_causality(section3_hpf, x, 500), ': section3_hpf might be causal')
print(submitted.test_causality(section4_rms, x, 500), ': section4_rms might be causal')
print(submitted.test_causality(section5_locsum, x[:,0], 500), ': section5_locsum might be causal')
print(submitted.test_causality(section6_threshold, x[:,0], 500), ': section6_threshold might be causal')
print(submitted.test_causality(section7_mindur, x[:,0], 500), ': section7_mindur might be causal')
print(submitted.test_causality(section8_lids, x[:,0], 500), ': section8_lids might be causal')

True : section1_lpf might be causal.
False : section3_hpf might be causal
True : section4_rms might be causal
False : section5_locsum might be causal
True : section6_threshold might be causal
True : section7_mindur might be causal
True : section8_lids might be causal

The results are:

• LPF, RMS, Threshold, MinDur and LIDS might be causal.
• HPF and Local Sum are non-causal.

Can you prove that RMS, Threshold, MinDur and LIDS are causal? Can you prove that HPF and Local Sum are non-causal? Can you prove using pencil and paper (not software) that Downsample is also causal?

6. Stability

A system is stable if every bounded input, $x[n]$, produces a bounded output $y[n]$. Bounded means that the maximum amplitude is finite. We usually write it something like this:

$$|x[n]|\le X_{max}<\infty \Rightarrow |y[n]|\le Y_{max}<\infty$$

No physical signal ever grows to infinity. What physical signals sometimes do, instead, is to continue growing without bound: The longer you wait, the larger the signal will grow.

Computer systems, on the other hand, eventually generate the largest number that the system can represent. Smart software like numpy will, at that point, generate a special symbol called np.inf that means "a number larger than I can represent." For example:

x = np.arange(100,1000,100)
y = np.exp(x)
for xn,yn in zip(x,y):
    print('exp(%d) ='%(xn), yn)
print('y[n] reaches a maximum absolute value of:',np.amax(np.abs(y)))

exp(100) = 2.6881171418161356e+43
exp(200) = 7.225973768125749e+86
exp(300) = 1.9424263952412558e+130
exp(400) = 5.221469689764144e+173
exp(500) = 1.4035922178528375e+217
exp(600) = 3.7730203009299397e+260
exp(700) = 1.0142320547350045e+304
exp(800) = inf
exp(900) = inf
y[n] reaches a maximum absolute value of: inf

/var/folders/9g/b0tx1gdj2vz6g_w7_h1w3wx80000gn/T/ipykernel_26845/9103472.py:2: RuntimeWarning: overflow encountered in exp
  y = np.exp(x)

Since numpy has this property, we can perform a simple test to see if any system is stable:

importlib.reload(submitted)
help(submitted.test_stability)

Help on function test_stability in module submitted:

test_stability(system, signal)
    Run a test, using the provided signal, to determine whether the system might be stable:
    If max(abs(signal)) is finite but max(abs(system(signal))) is infinite, then
    the system might be unstable, so this function should return False.
    If max(abs(signal)) is infinite or max(abs(system(signal))) is finite, then
    the system might be stable, so this function should return True.

• As with the other tests, a True result does not guarantee that the system is stable. It's possible that the system is unstable, but we just haven't found the signal that causes it to exhibit unstable behavior. If the test returns True with a long enough input signal, then it's plausible to think the system is stable, but it's not guaranteed.
• Unlike the other tests, a False result does not guarantee that the system is unstable! It's possible that the output signal is finite but extremely large -- so large that numpy can't represent it. In that case the system is stable, even though our test would erroneously call it unstable.

Revise "test_stability," and test it with the following block.

importlib.reload(submitted)

x = downsample_output

print(submitted.test_stability(section1_lpf, x), ': section1_lpf might be stable.')
print(submitted.test_stability(section2_downsample, x), ': section2_downsample might be stable')
print(submitted.test_stability(section3_hpf, x), ': section3_hpf might be stable')
print(submitted.test_stability(section4_rms, x), ': section4_rms might be stable')
print(submitted.test_stability(section5_locsum, x[:,0]), ': section5_locsum might be stable')
print(submitted.test_stability(section6_threshold, x[:,0]), ': section6_threshold might be stable')
print(submitted.test_stability(section7_mindur, x[:,0]), ': section7_mindur might be stable')
print(submitted.test_stability(section8_lids, x[:,0]), ': section8_lids might be stable')

True : section1_lpf might be stable.
True : section2_downsample might be stable
True : section3_hpf might be stable
True : section4_rms might be stable
True : section5_locsum might be stable
True : section6_threshold might be stable
True : section7_mindur might be stable
True : section8_lids might be stable

Well! It shouldn't surprise you very much to learn that all those systems are stable, since we ran them earlier without generating infinite outputs. Can you use pencil and paper to prove that any bounded input produces a bounded output from the Downsample, HPF, RMS, LocSum, Threshold, MinDur, and LIDS systems?

What happens if we run it with np.exp as the system?

print(submitted.test_stability(np.exp, x), ': np.exp might be stable if x is the input')
print(submitted.test_stability(np.exp, 10000000*x), ': np.exp might be stable if 10000000*x is the input')

True : np.exp might be stable if x is the input
False : np.exp might be stable if 10000000*x is the input

/Users/jhasegaw/Dropbox/mark/teaching/ece401/ece401labs/26fall/mp1/src/submitted.py:75: RuntimeWarning: overflow encountered in exp
  if np.amax(np.abs(signal)) == np.inf or np.amax(np.abs(system(signal))) != np.inf:

Well, that's interesting. Is np.exp a stable or unstable system? Can you use pencil and paper to figure it out?

Let's try another one:

def positive_feedback(x):
    y = np.zeros(x.shape)
    y[0] = x[0]
    for n in range(1,len(x)):
        y[n] = x[n] + 1.1 * y[n-1]
    return y

print(submitted.test_stability(positive_feedback, x), ': positive_feedback might be stable')

False : positive_feedback might be stable

/var/folders/9g/b0tx1gdj2vz6g_w7_h1w3wx80000gn/T/ipykernel_26845/1930980275.py:5: RuntimeWarning: overflow encountered in multiply
  y[n] = x[n] + 1.1 * y[n-1]

Stability is a somewhat counter-intuitive concept because it deals with infinities. The simple accumulator, $y[n]=x[n]+y[n-1]$, is unstable, but the output grows very slowly, so it would be hard to generate a long enough input signal to detect the instability. If we add the positive feedback to make $y[n]=x[n]+1.1y[n-1]$, however, the output grows much more quickly, so the instability is easily detected.

7. Submitting your answers to the autograder

Once your submitted.py gives correct answers to all of the blocks above, try running grade.py in the block below to see if your code passes. If it passes on your own computer, it will probably pass on the autograder; submit it online to find out. You can submit your code to the autograder as many times as you wish, until the deadline; the autograder will automatically keep your best score.

!pip install gradescope_utils

Requirement already satisfied: gradescope_utils in /opt/anaconda3/lib/python3.13/site-packages (0.5.0)

!python grade.py

........
----------------------------------------------------------------------
Ran 8 tests in 0.081s

OK