Process a low-rate measured energy signal
Inspect a daily-unit spectrum, an explicitly acausal filter and the limits of ten-minute sampling.
Goal and how to run
This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.
Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.
Setup and parameters
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from scipy.signal import welch, butter, sosfiltfilt, sosfreqz
from datetime import datetime
from IPython.display import display as _display
def display(value):
if isinstance(value,pd.DataFrame):
table=value.to_html(max_rows=12,escape=True)
style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
return _display(value)
RECIPE_ID='en03-energy-filtering'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
"axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6 Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads Core data mode: embedded snapshot; no network requests
Load the embedded snapshot
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
"""Read embedded data from this notebook, including a renamed copy."""
for path in sorted(Path.cwd().glob("*.ipynb")):
try:
notebook = json.loads(path.read_text(encoding="utf-8"))
metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
blob = metadata.get("snapshots", {}).get(snapshot_key)
if metadata.get("recipe_id") != recipe_id or not blob:
continue
if blob.get("sha256") != expected_sha256:
continue
raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
if hashlib.sha256(raw).hexdigest() != expected_sha256:
raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
return json.loads(raw)
except (OSError, json.JSONDecodeError):
continue
raise FileNotFoundError(
"Save/open the complete .ipynb in the active working folder. "
"The embedded snapshot could not be located for " + recipe_id)
SNAPSHOT_HASHES={'energy': '81117ee0a0706bf416bedc50f23b455ee8dd771e588c91b6a420dc3b184f638f'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['energy']
Read the method functions
def audit_numeric_table(frame,columns):
"""Audit observed rows; never silently impute or remove them."""
if frame.empty or frame.index.has_duplicates:
raise ValueError('Expected a nonempty table with unique row IDs')
missing=[name for name in columns if name not in frame.columns]
if missing:raise ValueError('Missing required columns: '+str(missing))
values=frame[columns].to_numpy(dtype=float)
if not np.isfinite(values).all():raise ValueError('Non-finite measurements; inspect missing values explicitly')
return {'rows':len(frame),'features':len(columns),'finite':True,'duplicate_measurement_rows':int(frame[columns].duplicated().sum())}
print("Method ready: audit_numeric_table")
Method ready: audit_numeric_table
def iso_datetimes(values):
"""Parse source ISO dates with the standard library; keep tables as ISO strings."""
return [datetime.fromisoformat(str(value)) for value in values]
print("Method ready: iso_datetimes")
Method ready: iso_datetimes
Data and model boundary
The building record contains real ten-minute energy and environmental observations from one source building. Airport weather fields include upstream interpolation and are not colocated calibration references. Beam, circuit and PID examples are declared numerical design models; computed responses do not establish hardware performance or commissioning acceptance.
Steps
1. Set the sampling scale in days
Ten-minute observations give 144 samples/day and Nyquist 72 cycles/day. They cannot resolve 50/60 Hz mains waveforms. This workflow analyzes variation of interval-average energy/power on much slower time scales.
energy=snapshot("energy");frame=pd.DataFrame(energy["rows"],columns=energy["columns"])
frame["timestamp"]=frame["date"].astype(str);frame=frame.sort_values("timestamp").reset_index(drop=True)
assert len(frame)==19735 and frame.timestamp.is_unique
dates=iso_datetimes(frame.timestamp)
assert all((b-a).total_seconds()==600 for a,b in zip(dates,dates[1:]))
audit=audit_numeric_table(frame,["Appliances","lights","T1","RH_1"])
assert frame.Appliances.ge(0).all() and frame.lights.ge(0).all()
frame["appliance_average_power_W"]=frame.Appliances/(600/3600)
print("Observed rows:",len(frame),"; energy: Wh/10 min; derived interval-average power: W; source timezone unspecified")
display(frame[["timestamp","Appliances","appliance_average_power_W","T1","RH_1"]].head(8))
power=frame.appliance_average_power_W.to_numpy(float);sampling_per_day=144.0;cutoff_cycles_per_day=6.0
sos=butter(4,cutoff_cycles_per_day,btype="lowpass",fs=sampling_per_day,output="sos")
smoothed=sosfiltfilt(sos,power)
fig,ax=plt.subplots(figsize=(10,4));count=7*144
ax.plot(iso_datetimes(frame.timestamp.iloc[:count]),power[:count],label="Observed interval-average power",lw=.7,alpha=.7)
ax.plot(iso_datetimes(frame.timestamp.iloc[:count]),smoothed[:count],label="Acausal low-pass processing",lw=1.4)
ax.set(xlabel="Source-local time",ylabel="Average / processed power (W)",title="Measured source and offline processing");ax.legend();plt.tight_layout();plt.show()
Observed rows: 19735 ; energy: Wh/10 min; derived interval-average power: W; source timezone unspecified
| timestamp | Appliances | appliance_average_power_W | T1 | RH_1 | |
|---|---|---|---|---|---|
| 0 | 2016-01-11 17:00:00 | 60 | 360.0 | 19.890000 | 47.596667 |
| 1 | 2016-01-11 17:10:00 | 60 | 360.0 | 19.890000 | 46.693333 |
| 2 | 2016-01-11 17:20:00 | 50 | 300.0 | 19.890000 | 46.300000 |
| 3 | 2016-01-11 17:30:00 | 50 | 300.0 | 19.890000 | 46.066667 |
| 4 | 2016-01-11 17:40:00 | 60 | 360.0 | 19.890000 | 46.333333 |
| 5 | 2016-01-11 17:50:00 | 50 | 300.0 | 19.890000 | 46.026667 |
| 6 | 2016-01-11 18:00:00 | 60 | 360.0 | 19.890000 | 45.766667 |
| 7 | 2016-01-11 18:10:00 | 60 | 360.0 | 19.856667 | 45.560000 |
2. Estimate the low-rate spectrum
With fs expressed as samples/day, Welch frequencies are cycles/day and density is W² per (cycles/day). A daily-pattern peak is a description of this source period, not evidence of high-rate electrical harmonics.
frequency,psd=welch(power-power.mean(),fs=sampling_per_day,nperseg=7*144,noverlap=3*144)
positive=frequency>0
fig,ax=plt.subplots(figsize=(8,5));ax.loglog(frequency[positive],psd[positive])
ax.axvline(cutoff_cycles_per_day,color="red",linestyle="--",label="Declared cutoff")
ax.set(xlabel="Frequency (cycles/day)",ylabel="Power PSD (W² per cycles/day)",title="Welch density at the actual observation cadence");ax.legend();plt.tight_layout();plt.show()
3. Read the double-pass response
Forward/backward processing is acausal and squares the single-pass magnitude. It is useful for offline exploration, but cannot be deployed unchanged as a real-time controller filter. Boundaries and ringing can affect the processed series.
response_frequency,response=sosfreqz(sos,worN=2048,fs=sampling_per_day)
fig,ax=plt.subplots(figsize=(8,4));ax.plot(response_frequency,abs(response),label="Single pass")
ax.plot(response_frequency,abs(response)**2,label="Forward/backward magnitude")
ax.axvline(cutoff_cycles_per_day,color="gray",linestyle="--")
ax.set(xlabel="Frequency (cycles/day)",ylabel="Magnitude ratio",xlim=(0,20),title="Declared offline filter response");ax.legend();plt.tight_layout();plt.show()
4. Check FFT normalization and preserve original data
Use the full complex FFT for Parseval, avoiding an even/odd real-FFT endpoint mistake on 19735 observations. Keep measured and processed series separate and report any negative processed values without hiding them.
centered=power-power.mean();transform=np.fft.fft(centered)
assert np.isclose(np.sum(abs(transform)**2)/len(power)**2,np.mean(centered**2),rtol=1e-12)
assert np.isfinite(smoothed).all() and np.all(psd>=0)
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"timestamp":frame.timestamp,"observed_interval_average_W":power,"acausal_lowpass_W":smoothed}).to_csv(OUTPUT_DIR/"power-processing.csv",index=False)
pd.DataFrame({"cycles_per_day":frequency,"PSD_W_squared_per_cycles_per_day":psd}).to_csv(OUTPUT_DIR/"energy-spectrum.csv",index=False)
print("Cadence/FFT checks passed; negative processed values:",int(np.sum(smoothed<0)))
Cadence/FFT checks passed; negative processed values: 4
Exercises and reference explanation
- What is the Nyquist limit in cycles/day?
- Why square the single-pass magnitude?
Reference explanation. It is 72 cycles/day at 144 samples/day. Forward/backward filtering composes two passes and is acausal; its magnitude is the square of the one-pass magnitude.
Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.
Optional AI coding exercises
Prompt 1
Using NumPy/SciPy/Matplotlib, compare 3 and 6 cycles/day cutoffs without changing cadence labels.
Prompt 2
With the same packages, plot first/last-week boundary behavior and preserve processed versus measured values.
Human acceptance. FFT normalization passes for the odd observation count, spectra have stated units and filtering is labelled offline/acausal.
Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.
Optional online extension
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://archive.ics.uci.edu/dataset/374/appliances+energy+prediction',timeout=20)
response.raise_for_status()
print("Source:",response.url)
print(response.text[:700])
except requests.RequestException as error:
print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged
Sources, snapshots and next steps
- UCI energy — CC BY 4.0; credit UCI and original dataset authors; retrieved 2026-10-09. SHA-256
2fccf354445d886e7917620b0195db1f3e3e34d5a067a93b844694a4c561255a. Extract energydata_complete.csv from original download; preserve all rows and declared units; no imputation. Candanedo, L. (2017). Appliances Energy Prediction. UCI Machine Learning Repository. DOI: 10.24432/C5VC8G.
Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.