1. 項(xiàng)目概述什么是“最大8小時內(nèi)滑動平均”在數(shù)據(jù)分析、環(huán)境監(jiān)測、金融量化等領(lǐng)域我們常常需要處理時間序列數(shù)據(jù)。一個典型的需求是從連續(xù)的時間序列中找出任意連續(xù)的8小時窗口內(nèi)某個指標(biāo)比如PM2.5濃度、股票成交量、服務(wù)器負(fù)載的平均值然后從所有可能的8小時窗口中找出那個最大的平均值。這就是“計(jì)算最大8小時內(nèi)滑動平均”的核心任務(wù)。聽起來簡單但手動計(jì)算幾乎不可能尤其是面對動輒數(shù)萬、數(shù)十萬條記錄的數(shù)據(jù)時。MATLAB作為強(qiáng)大的數(shù)值計(jì)算和工程仿真平臺其向量化操作和豐富的內(nèi)置函數(shù)讓這類計(jì)算變得高效而優(yōu)雅。今天我就以一個從業(yè)多年的數(shù)據(jù)分析師視角帶你從零開始手把手實(shí)現(xiàn)這個功能并深入探討其中的技術(shù)細(xì)節(jié)、性能優(yōu)化和那些“教科書上不會寫”的避坑經(jīng)驗(yàn)。2. 核心思路與方案選型為什么是“滑動窗口”在動手寫代碼之前我們必須先理清思路。計(jì)算“最大8小時滑動平均”本質(zhì)上是一個滑動窗口統(tǒng)計(jì)問題。這里的“8小時”是窗口的寬度而“滑動”意味著這個窗口會沿著時間軸以一個固定的步長通常是1個數(shù)據(jù)點(diǎn)即逐點(diǎn)滑動移動每移動一次就計(jì)算一次窗口內(nèi)數(shù)據(jù)的平均值。2.1 方案對比循環(huán) vs. 向量化面對這個問題新手最容易想到的方法是使用for循環(huán)遍歷時間序列的每一個可能起點(diǎn)截取接下來8小時的數(shù)據(jù)計(jì)算平均值然后更新最大值。這個方法直觀但效率是硬傷。MATLAB的for循環(huán)在處理大規(guī)模數(shù)據(jù)時性能較差尤其是在腳本中直接操作時。向量化操作是MATLAB的靈魂。我們的目標(biāo)是盡可能利用MATLAB內(nèi)置的、用C/C優(yōu)化過的函數(shù)避免顯式循環(huán)。對于滑動平均MATLAB提供了幾個潛在的“武器”movmean函數(shù)這是最直接的工具專門用于計(jì)算移動平均值。語法簡潔性能優(yōu)異。卷積操作利用conv函數(shù)與一個全為1的向量進(jìn)行卷積再除以窗口長度可以實(shí)現(xiàn)滑動求和進(jìn)而得到平均值。這是一種更底層、更靈活的方法。filter函數(shù)作為信號處理工具箱的一員它本質(zhì)上也是在實(shí)現(xiàn)卷積可以用于計(jì)算滑動平均。對于“最大8小時滑動平均”這個具體任務(wù)movmean函數(shù)是首選。因?yàn)樗Z義清晰無需自己處理邊界條件如窗口在數(shù)據(jù)開頭和結(jié)尾時數(shù)據(jù)不足的問題并且經(jīng)過了高度優(yōu)化。2.2 數(shù)據(jù)準(zhǔn)備與關(guān)鍵假設(shè)在編碼前我們必須明確幾個前提這直接關(guān)系到代碼的健壯性時間間隔你的數(shù)據(jù)必須是等時間間隔的。例如每小時一個數(shù)據(jù)點(diǎn)或者每分鐘一個數(shù)據(jù)點(diǎn)。如果數(shù)據(jù)間隔不均勻直接滑動平均沒有物理意義需要先進(jìn)行重采樣或插值。窗口單位的轉(zhuǎn)換“8小時”是一個時間長度而你的數(shù)據(jù)索引通常是數(shù)據(jù)點(diǎn)序號。你需要知道數(shù)據(jù)的采樣頻率。例如如果你的數(shù)據(jù)是每小時一個點(diǎn)那么8小時窗口就對應(yīng)8個數(shù)據(jù)點(diǎn)。如果是每5分鐘一個點(diǎn)那么8小時窗口就對應(yīng)8 * 60 / 5 96個數(shù)據(jù)點(diǎn)。邊界處理movmean函數(shù)默認(rèn)會處理邊界。對于窗口起始部分不足8個點(diǎn)的情況它會計(jì)算已有數(shù)據(jù)的平均值這被稱為‘shrink’模式。我們需要決定這個行為是否符合需求。在環(huán)境標(biāo)準(zhǔn)計(jì)算中如計(jì)算“日最大8小時平均”通常要求窗口必須完整包含8個數(shù)據(jù)點(diǎn)不足的則不予計(jì)算。movmean可以通過指定‘Endpoints’參數(shù)來控制。注意如果你的數(shù)據(jù)時間戳是datetime格式而數(shù)值是單獨(dú)數(shù)組你需要先將時間戳轉(zhuǎn)換為等間隔的索引或者利用時間戳直接邏輯索引來構(gòu)造窗口。本文假設(shè)你已經(jīng)有了一個等間隔的數(shù)值向量data和對應(yīng)的采樣頻率Fs單位點(diǎn)/小時。3. 核心實(shí)現(xiàn)與代碼逐行解析理論清晰后我們進(jìn)入實(shí)戰(zhàn)環(huán)節(jié)。我將提供兩個版本的代碼一個基礎(chǔ)通用版一個考慮環(huán)境監(jiān)測實(shí)際應(yīng)用的加強(qiáng)版。3.1 基礎(chǔ)通用版實(shí)現(xiàn)假設(shè)我們有一個向量data它是按小時采樣的濃度數(shù)據(jù)例如PM2.5單位μg/m3。我們要找出任意連續(xù)8小時內(nèi)的最大平均濃度。% 基礎(chǔ)版計(jì)算最大8小時滑動平均 % 假設(shè)數(shù)據(jù) data 是每小時一個點(diǎn)的濃度序列 Fs 1; % 采樣頻率1點(diǎn)/小時 windowLengthHours 8; windowLengthPoints windowLengthHours * Fs; % 窗口長度 8個點(diǎn) % 使用 movmean 計(jì)算滑動平均 % ‘Endpoints’, ‘discard’ 表示在數(shù)據(jù)兩端當(dāng)窗口不能完整覆蓋時結(jié)果中丟棄這些位置的值。 % 這符合“必須完整8小時”的常見要求。 movingAvg movmean(data, windowLengthPoints, ‘Endpoints’, ‘discard’); % 找出所有滑動平均值中的最大值 max8hrAvg max(movingAvg); % 可選找出最大值發(fā)生的位置窗口的起始索引 [maxValue, maxIndex] max(movingAvg); % 注意maxIndex 對應(yīng)的是 movingAvg 向量中的位置。 % 要找到原始數(shù)據(jù)中對應(yīng)窗口的起始索引因?yàn)椤甦iscard’了前(windowLengthPoints-1)個點(diǎn)所以需要加上偏移量。 windowStartIndex maxIndex; % 因?yàn)閬G棄了端點(diǎn)所以 movingAvg 的第一個值對應(yīng)原始數(shù)據(jù)中第一個完整窗口的開始 fprintf(‘最大8小時滑動平均值為%.2f\n’, max8hrAvg); fprintf(‘該最大值出現(xiàn)在從第%d小時開始的8小時窗口內(nèi)。\n’, windowStartIndex);代碼解析與注意事項(xiàng)movmean參數(shù)詳解data: 輸入的時間序列向量。windowLengthPoints: 窗口長度以數(shù)據(jù)點(diǎn)數(shù)為單位。這里是8。‘Endpoints’, ‘discard’: 這是關(guān)鍵參數(shù)。它指定了在數(shù)據(jù)序列的開始和結(jié)尾當(dāng)滑動窗口無法被數(shù)據(jù)完全填滿時如何處理輸出。‘discard’會直接忽略這些不完整的窗口不在movingAvg中輸出它們的值。這對于尋找“完整8小時窗口內(nèi)的最大平均”至關(guān)重要。如果使用默認(rèn)值或‘shrink’則會用已有數(shù)據(jù)計(jì)算平均值可能導(dǎo)致結(jié)果偏大或偏小。索引對齊的坑 計(jì)算出的maxIndex是movingAvg向量中的索引。由于我們丟棄了前7個不完整窗口對于8點(diǎn)窗口movingAvg(1)實(shí)際上對應(yīng)的是原始數(shù)據(jù)data(1:8)這個完整窗口的平均值。因此原始數(shù)據(jù)中對應(yīng)窗口的起始索引就是maxIndex。如果你使用了‘shrink’或其他端點(diǎn)處理方法這個對應(yīng)關(guān)系會發(fā)生變化必須仔細(xì)推算。3.2 環(huán)境監(jiān)測應(yīng)用加強(qiáng)版在實(shí)際環(huán)境空氣質(zhì)量評價(jià)中“日最大8小時平均”是一個重要指標(biāo)。它的規(guī)則更具體一天有24小時但計(jì)算的是“移動的8小時平均”即從0點(diǎn)到24點(diǎn)每一個小時作為起始點(diǎn)取其后8小時的平均值全天共有17個24-81這樣的滑動平均值再取這17個值中的最大值作為當(dāng)日的“日最大8小時平均”。并且通常要求每小時的數(shù)據(jù)是有效的非缺失。下面我們模擬一個更真實(shí)的場景數(shù)據(jù)包含日期時間信息并且可能存在缺失值用NaN表示。% 加強(qiáng)版考慮日期時間和缺失值計(jì)算“日最大8小時平均” % 1. 生成模擬數(shù)據(jù)假設(shè)為2023年某一天每小時的數(shù)據(jù) dateVector datetime(2023, 6, 1, 0, 0, 0):hours(1):datetime(2023, 6, 1, 23, 0, 0); % 模擬一些隨機(jī)濃度數(shù)據(jù)并插入一些缺失值(NaN) rng(‘default’); % 保證可重復(fù)性 data 30 20 * randn(size(dateVector)); % 均值為50的正態(tài)分布隨機(jī)數(shù) data([5, 15, 22]) NaN; % 在第5, 15, 22小時設(shè)置數(shù)據(jù)缺失 % 2. 處理缺失值 - 對于滑動平均常見的簡單處理是線性插值 dataFilled fillmissing(data, ‘linear’); % 使用線性插值填充NaN % 注意也可以使用 ‘previous’, ‘next’ 或 ‘nearest’。選擇取決于實(shí)際業(yè)務(wù)邏輯。 % 如果缺失值過多插值可能引入較大誤差需要評估。 % 3. 計(jì)算8小時滑動平均完整窗口 windowHours 8; movingAvgFull movmean(dataFilled, windowHours, ‘Endpoints’, ‘discard’); % 4. 找出最大值及其位置 [maxAvgValue, maxAvgIdxInMoving] max(movingAvgFull); % 計(jì)算該最大值對應(yīng)的原始數(shù)據(jù)時間窗口的起始時間 % movingAvgFull(1) 對應(yīng)原始時間 dateVector(1) 到 dateVector(8) 的平均值 windowStartTime dateVector(maxAvgIdxInMoving); windowEndTime dateVector(maxAvgIdxInMoving windowHours - 1); % 結(jié)束時間是起始時間7小時 % 5. 輸出結(jié)果 fprintf(‘日期%s\n’, datestr(dateVector(1), ‘yyyy-mm-dd’)); fprintf(‘經(jīng)過線性插值處理后日最大8小時平均濃度為%.2f μg/m3\n’, maxAvgValue); fprintf(‘該最大值對應(yīng)的8小時窗口為%s 至 %s\n’, ... datestr(windowStartTime, ‘HH:MM’), datestr(windowEndTime, ‘HH:MM’)); % 6. 可視化繪制原始數(shù)據(jù)、插值后數(shù)據(jù)及滑動平均曲線 figure(‘Position’, [100, 100, 1200, 500]); subplot(2,1,1); plot(dateVector, data, ‘o-‘, ‘DisplayName’, ‘原始數(shù)據(jù)含NaN’); hold on; plot(dateVector, dataFilled, ‘x–‘, ‘DisplayName’, ‘插值后數(shù)據(jù)’); xlabel(‘時間’); ylabel(‘濃度 (μg/m3)’); title(‘原始數(shù)據(jù)與缺失值處理’); legend(‘Location’, ‘best’); grid on; subplot(2,1,2); % 為滑動平均結(jié)果生成對應(yīng)的時間軸丟棄了前7個點(diǎn) timeForMovingAvg dateVector(1:end-windowHours1) hours((windowHours-1)/2); % 將時間點(diǎn)標(biāo)在窗口中部 plot(timeForMovingAvg, movingAvgFull, ‘s-‘, ‘LineWidth’, 1.5, ‘DisplayName’, ‘8小時滑動平均’); hold on; % 標(biāo)記出最大值點(diǎn) plot(timeForMovingAvg(maxAvgIdxInMoving), maxAvgValue, ‘r*’, ‘MarkerSize’, 15, ‘DisplayName’, ‘日最大8小時平均’); xlabel(‘時間窗口中心點(diǎn)’); ylabel(‘平均濃度 (μg/m3)’); title(‘8小時滑動平均序列與最大值’); legend(‘Location’, ‘best’); grid on;關(guān)鍵點(diǎn)解析與實(shí)操心得缺失值處理是重中之重movmean函數(shù)遇到NaN時整個窗口的平均值也會是NaN。這會導(dǎo)致最大值查找失敗max函數(shù)會忽略NaN但你可能得到的是一個非完整窗口的最大值。因此必須先處理缺失值。fillmissing函數(shù)非常強(qiáng)大‘linear’插值適用于連續(xù)變化的數(shù)據(jù)。但在實(shí)際業(yè)務(wù)中需要根據(jù)數(shù)據(jù)缺失機(jī)制和行業(yè)規(guī)范選擇方法有時甚至需要將缺失過多的小時所在日的計(jì)算視為無效。時間戳對齊滑動平均結(jié)果movingAvgFull的長度比原始數(shù)據(jù)短。為了繪圖或分析需要為其創(chuàng)建正確的時間標(biāo)簽。常見的做法是將平均值對應(yīng)的時間點(diǎn)放在窗口的中間時刻如上例代碼所示這樣在圖上看起來更合理。而查找出的maxAvgIdxInMoving對應(yīng)的是這個“中間時刻”序列的索引要反推回窗口的起止時間需要做簡單的加減運(yùn)算。‘Endpoints’, ‘discard’的必然性在環(huán)境標(biāo)準(zhǔn)計(jì)算中必須使用此參數(shù)。因?yàn)橐惶靸啥说拇翱谌?-7點(diǎn)17-24點(diǎn)是不完整的24小時內(nèi)的8小時窗口不符合“日內(nèi)滑動”的定義。計(jì)算時只考慮從0點(diǎn)至16點(diǎn)開始的共17個完整窗口。4. 性能優(yōu)化與高級技巧當(dāng)數(shù)據(jù)量極大例如多年、多站點(diǎn)的每小時數(shù)據(jù)時基礎(chǔ)方法可能仍有優(yōu)化空間。此外一些特殊需求也需要更靈活的方案。4.1 處理超長序列與分塊計(jì)算對于長達(dá)數(shù)年的每小時數(shù)據(jù)直接計(jì)算內(nèi)存占用可能很高。雖然movmean已經(jīng)優(yōu)化得很好但我們可以考慮分日計(jì)算因?yàn)椤叭兆畲?小時平均”本身就是按日統(tǒng)計(jì)的。% 假設(shè)我們有長時間序列數(shù)據(jù) dates 和 values % 首先將數(shù)據(jù)按日期分組 [year, month, day] ymd(dates); % 需要 datetime 數(shù)組 dateGroups findgroups(year, month, day); % 為每一天創(chuàng)建一個分組ID % 預(yù)分配結(jié)果數(shù)組 uniqueDates unique(dates, ‘day’); % 獲取不重復(fù)的日期 max8hrDaily zeros(size(uniqueDates)); % 對每一天進(jìn)行循環(huán)計(jì)算 for i 1:length(uniqueDates) dayMask dates uniqueDates(i) dates uniqueDates(i) days(1); dataOfDay values(dayMask); % 處理缺失值這里簡單用前后值均值填充實(shí)際需謹(jǐn)慎 dataFilled fillmissing(dataOfDay, ‘linear’); % 確保一天有24個數(shù)據(jù)點(diǎn)處理可能的嚴(yán)重缺失 if length(dataFilled) 24 movingAvg movmean(dataFilled, 8, ‘Endpoints’, ‘discard’); max8hrDaily(i) max(movingAvg); else max8hrDaily(i) NaN; % 數(shù)據(jù)不全記為缺失 end end % 現(xiàn)在 max8hrDaily 就是每一天的“日最大8小時平均”4.2 自定義滑動窗口函數(shù)以應(yīng)對復(fù)雜邏輯如果業(yè)務(wù)邏輯非常特殊比如窗口長度可變或者計(jì)算的不是算術(shù)平均而是其他統(tǒng)計(jì)量如中位數(shù)、百分位數(shù)可以自定義滑動窗口函數(shù)。% 示例計(jì)算8小時滑動中位數(shù)對異常值更魯棒 windowLen 8; data randn(1000,1); % 模擬數(shù)據(jù) % 方法使用循環(huán)但利用預(yù)分配和向量索引提高效率 n length(data); result zeros(n - windowLen 1, 1); % 預(yù)分配結(jié)果數(shù)組 for startIdx 1:(n - windowLen 1) windowData data(startIdx : startIdx windowLen - 1); result(startIdx) median(windowData); end maxSlidingMedian max(result);雖然用了循環(huán)但對于窗口操作MATLAB R2016a以后版本對for循環(huán)進(jìn)行了JIT即時編譯加速在不是極端性能瓶頸的場景下這種寫法清晰易懂。當(dāng)然也可以探索用arrayfun或編寫MEX文件來進(jìn)一步優(yōu)化。4.3 利用卷積conv實(shí)現(xiàn)底層滑動平均理解movmean的底層原理有助于解決更復(fù)雜的問題。滑動平均可以通過卷積實(shí)現(xiàn)windowLen 8; kernel ones(windowLen, 1) / windowLen; % 卷積核長度為8每個元素為1/8 movingAvgConv conv(data, kernel, ‘valid’); % ‘valid’模式只返回完全重疊的部分相當(dāng)于‘discard’‘valid’模式的結(jié)果長度是length(data) - windowLen 1與movmean(data, windowLen, ‘Endpoints’, ‘discard’)結(jié)果完全相同。這種方法在你想自定義加權(quán)平均如指數(shù)加權(quán)時特別有用只需修改kernel向量即可。5. 常見問題、錯誤排查與調(diào)試技巧在實(shí)際操作中你幾乎一定會遇到下面這些問題。這里是我的“踩坑”實(shí)錄和解決方案。5.1 數(shù)據(jù)長度與窗口長度不匹配問題計(jì)算時MATLAB報(bào)錯“窗口長度必須小于或等于輸入長度”或結(jié)果的長度出乎意料。排查檢查你的windowLengthPoints計(jì)算是否正確。確保它是標(biāo)量整數(shù)。檢查輸入數(shù)據(jù)data是否是向量。movmean也支持矩陣按指定維度計(jì)算如果data是矩陣需要指定維度參數(shù)如movmean(data, k, 1)對列滑動。回想‘Endpoints’參數(shù)的影響。如果使用‘discard’輸出長度會是length(data) - windowLengthPoints 1。如果你期望輸出長度與輸入相同應(yīng)使用‘shrink’默認(rèn)或‘fill’。5.2 結(jié)果全是NaN或包含NaN問題計(jì)算出的movingAvg里有很多甚至全部是NaN。原因與解決輸入數(shù)據(jù)包含NaN這是最常見原因。使用any(isnan(data))檢查。務(wù)必在計(jì)算前處理缺失值插值、刪除或標(biāo)記。窗口內(nèi)全是NaN即使做了插值如果數(shù)據(jù)開頭或結(jié)尾連續(xù)缺失插值可能失敗。考慮使用‘omitnan’選項(xiàng)MATLAB R2015b以上movmean(data, k, ‘omitnan’)。這個選項(xiàng)會在計(jì)算每個窗口的平均值時忽略該窗口內(nèi)的NaN。但要注意這可能導(dǎo)致窗口實(shí)際用于計(jì)算的數(shù)據(jù)點(diǎn)數(shù)少于k從而影響結(jié)果的可比性需結(jié)合業(yè)務(wù)判斷。5.3 最大值對應(yīng)的時間窗口找錯問題找到了最大平均值但根據(jù)索引回溯到原始數(shù)據(jù)時發(fā)現(xiàn)對應(yīng)的8小時窗口不對。調(diào)試步驟打印關(guān)鍵索引在計(jì)算后立即打印maxIndex,length(data),length(movingAvg)確認(rèn)它們的關(guān)系。手動驗(yàn)證一個小例子用一個人工構(gòu)造的簡單數(shù)組如data [1:24]運(yùn)行你的代碼。因?yàn)榈炔顢?shù)列的平均值就是中間值你可以很容易地心算出最大8小時平均應(yīng)該是[17:24]這個窗口平均值為20.5。檢查你的程序結(jié)果是否匹配。繪制示意圖像前面的加強(qiáng)版代碼一樣將原始數(shù)據(jù)、滑動平均序列以及標(biāo)記的最大值點(diǎn)畫在同一張圖上。視覺檢查是最有效的調(diào)試手段之一。5.4 處理非整點(diǎn)或不等間隔數(shù)據(jù)問題數(shù)據(jù)時間戳不是規(guī)整的整點(diǎn)時間或者采樣間隔不穩(wěn)定。解決方案重采樣使用retime針對timetable或resample針對信號函數(shù)將數(shù)據(jù)插值或聚合到等間隔的時間網(wǎng)格上例如每小時一個點(diǎn)。這是最規(guī)范的做法。基于時間戳的滑動如果堅(jiān)持使用原始時間戳你需要編寫自定義循環(huán)。在循環(huán)中對于每一個數(shù)據(jù)點(diǎn)i使用邏輯索引找出時間在[time(i), time(i)hours(8)]范圍內(nèi)的所有數(shù)據(jù)點(diǎn)然后計(jì)算它們的平均值。這種方法計(jì)算量很大但能最大程度保留原始信息。% 偽代碼示意 times ... % datetime 向量 values ... % 數(shù)值向量 maxAvg -inf; for i 1:length(times) windowMask times times(i) times times(i) hours(8); if sum(windowMask) 6 % 至少需要一定數(shù)量的數(shù)據(jù)點(diǎn)例如6個 avg mean(values(windowMask), ‘omitnan’); if avg maxAvg maxAvg avg; bestStartTime times(i); end end end5.5 內(nèi)存不足Out of Memory問題處理超大型數(shù)組時MATLAB報(bào)內(nèi)存錯誤。優(yōu)化策略使用單精度如果數(shù)據(jù)精度要求不高在數(shù)據(jù)導(dǎo)入時使用single類型data single(yourData);可以減半內(nèi)存占用。分塊處理如4.1節(jié)所示將數(shù)據(jù)按天、按月分割處理每次只加載一部分到內(nèi)存。避免創(chuàng)建中間大數(shù)組例如movmean的結(jié)果是一個新數(shù)組。如果原始數(shù)據(jù)很大這個結(jié)果數(shù)組也很大。如果后續(xù)只需要最大值可以考慮分塊計(jì)算并實(shí)時比較更新最大值而不是保存整個滑動平均序列。使用內(nèi)存映射文件對于存儲在磁盤上的巨型數(shù)據(jù)文件可以使用memmapfile函數(shù)進(jìn)行內(nèi)存映射實(shí)現(xiàn)按需訪問而不是一次性全部讀入。6. 擴(kuò)展應(yīng)用從滑動平均到滑動統(tǒng)計(jì)掌握了滑動平均你就可以輕松擴(kuò)展到其他滑動窗口統(tǒng)計(jì)量MATLAB提供了統(tǒng)一的movXXX函數(shù)家族movmedian: 滑動中位數(shù)抗噪聲movstd: 滑動標(biāo)準(zhǔn)差看波動movvar: 滑動方差movsum: 滑動總和movmin/movmax: 滑動最小/最大值例如在金融分析中我們常看股價(jià)的20日滑動標(biāo)準(zhǔn)差波動率在工業(yè)監(jiān)控中看設(shè)備溫度最近1小時的滑動最大值是否超閾值。其調(diào)用語法與movmean高度一致。最后關(guān)于工具版本我強(qiáng)烈建議使用MATLAB R2016a或更高版本。這些版本對movmean等函數(shù)以及循環(huán)的JIT編譯都有了顯著優(yōu)化性能提升非常明顯。如果你還在使用更舊的版本升級帶來的效率提升可能會讓你驚喜。