处理低采样率实测能耗信号

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

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']

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

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

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()

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()

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

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