1. 項目緣起為什么需要計算8小時滑動平均在數(shù)據(jù)分析、信號處理和環(huán)境監(jiān)測等領(lǐng)域我們常常會遇到一種需求評估某個指標(biāo)在特定時間窗口內(nèi)的平均表現(xiàn)并且這個窗口需要像“滑尺”一樣隨著時間推移而移動。滑動平均或者說移動平均就是處理這類需求的經(jīng)典工具。它能夠有效平滑數(shù)據(jù)中的短期波動和噪聲幫助我們更清晰地觀察數(shù)據(jù)的長期趨勢和周期性變化。那么為什么是“8小時”這個窗口這個需求在實際工作中非常普遍。一個典型的場景是空氣質(zhì)量評價。許多國家和地區(qū)的環(huán)境標(biāo)準(zhǔn)例如對臭氧、PM2.5等污染物的評價不僅看日均濃度更會關(guān)注“日最大8小時平均濃度”。這個指標(biāo)的計算方法是以一天中每一個小時作為結(jié)束點向前追溯8個小時計算這8個小時的濃度平均值然后從全天24個這樣的8小時平均值中找出最大的那一個。這個最大值才是評價當(dāng)天空氣質(zhì)量是否超標(biāo)的關(guān)鍵依據(jù)。類似地在工業(yè)過程控制、金融數(shù)據(jù)分析如計算某只股票在過去N個交易日的平均價格中滑動平均也是基礎(chǔ)但核心的操作。手動計算這個值非常繁瑣尤其是處理長時間序列數(shù)據(jù)時。MATLAB作為強(qiáng)大的數(shù)值計算和數(shù)據(jù)分析環(huán)境自然是完成這項任務(wù)的利器。但MATLAB本身并沒有一個直接叫做“max_8hr_moving_average”的函數(shù)。我們需要利用其強(qiáng)大的數(shù)組操作和函數(shù)組合能力來構(gòu)建一個高效、可靠的計算流程。這不僅僅是調(diào)用一個函數(shù)那么簡單它涉及到對滑動窗口概念的深刻理解、對MATLAB向量化編程的熟練運用以及對邊界條件、計算效率等細(xì)節(jié)的考量。2. 滑動平均的核心原理與MATLAB實現(xiàn)思路在動手寫代碼之前我們必須先吃透滑動平均的數(shù)學(xué)本質(zhì)。對于一個離散的時間序列數(shù)據(jù)y [y1, y2, y3, ..., yN]窗口長度為w在本文中w8的滑動平均會生成一個新的序列MA。對于新序列中的第i個元素MA(i)其計算公式為MA(i) (y(i-w1) y(i-w2) ... y(i)) / w其中i的取值范圍是從w到N。也就是說第一個有效的滑動平均值是從原始數(shù)據(jù)的第w個數(shù)據(jù)點開始計算的它代表了前w個數(shù)據(jù)的平均值。這里就引出了兩個關(guān)鍵問題邊界處理對于序列開頭i w的部分沒有足夠的數(shù)據(jù)填滿窗口這些位置的傳統(tǒng)滑動平均值是未定義的。在MATLAB中我們常見的處理方式是返回NaN非數(shù)字或者使用較小的窗口如1到i進(jìn)行計算。在環(huán)境標(biāo)準(zhǔn)計算中通常要求嚴(yán)格滿足8小時窗口因此前7個位置對于小時數(shù)據(jù)的結(jié)果應(yīng)為NaN或被視為無效。計算效率最直觀的方法是寫一個循環(huán)為每個i計算一次窗口內(nèi)數(shù)據(jù)的和。當(dāng)數(shù)據(jù)量N很大時例如多年的小時數(shù)據(jù)這種方法的效率很低。我們需要利用MATLAB的向量化操作來提升性能?；谝陨显鞰ATLAB中有幾種主流的實現(xiàn)思路思路一使用movmean函數(shù)最簡單直接這是R2016a版本后引入的官方函數(shù)專為移動計算設(shè)計。其基本語法是M movmean(A, k)其中k是窗口長度。對于我們的需求可以寫作window_size 8; eight_hr_ma movmean(hourly_data, [window_size-1, 0], ‘Endpoints’, ‘fill’);這里[window_size-1, 0]指定了一個非對稱窗口包含當(dāng)前點之前的7個點和當(dāng)前點本身總共8個點。‘Endpoints’, ‘fill’參數(shù)指定在數(shù)據(jù)開頭不足窗口長度的地方用NaN填充。這完美契合了環(huán)境標(biāo)準(zhǔn)中“向前追溯8小時”的定義。得到8小時滑動平均序列后再用max函數(shù)忽略NaN找出最大值即可。思路二使用卷積操作conv理解本質(zhì)靈活性強(qiáng)卷積是信號處理中實現(xiàn)滑動平均的數(shù)學(xué)基礎(chǔ)。一個長度為w的滑動平均等價于與一個元素全為1/w的長度為w的向量進(jìn)行卷積。window_size 8; kernel ones(window_size, 1) / window_size; % 創(chuàng)建平均核 eight_hr_ma_conv conv(hourly_data, kernel, ‘valid’);使用‘valid’模式時conv只計算那些不需要補零的部分其結(jié)果長度是N - window_size 1。它直接從第8個點開始輸出有效平均值前7個點被“丟棄”了。這同樣符合標(biāo)準(zhǔn)但需要注意結(jié)果序列與原始序列索引的對應(yīng)關(guān)系。思路三使用循環(huán)與向量化優(yōu)化深入控制教學(xué)意義為了深入理解過程我們可以從循環(huán)開始然后優(yōu)化。最樸素的循環(huán)寫法N length(hourly_data); window_size 8; eight_hr_ma_loop zeros(N, 1) * NaN; % 預(yù)先填充NaN for i window_size:N window_data hourly_data(i-window_size1:i); eight_hr_ma_loop(i) mean(window_data); end這個循環(huán)清晰易懂但速度慢。一個經(jīng)典的向量化優(yōu)化是使用累積和。我們先計算原始序列的累積和cumsum那么任意區(qū)間[i, j]的和就可以通過cumsum(j) - cumsum(i-1)快速得到。cs [0; cumsum(hourly_data)]; % 在開頭補一個0方便計算 eight_hr_ma_cumsum (cs(window_size1:end) - cs(1:end-window_size)) / window_size; % 結(jié)果長度為 N - window_size 1需要前面補上 NaN 以對齊原始序列 eight_hr_ma_final [nan(window_size-1, 1); eight_hr_ma_cumsum];這種方法在數(shù)學(xué)上等價于卷積但通過累積和避免了卷積函數(shù)可能的一些額外開銷在處理超大型數(shù)據(jù)時有時性能更優(yōu)。注意在計算“日最大8小時平均”時原始數(shù)據(jù)通常是按小時排列的濃度值。你需要確保數(shù)據(jù)是連續(xù)的沒有缺失的小時。如果存在缺失上述方法都會將其視為一個有效數(shù)據(jù)可能是0或NaN參與計算導(dǎo)致結(jié)果錯誤。因此數(shù)據(jù)預(yù)處理如插值或標(biāo)記缺失是必不可少的先決步驟。3. 從原理到實踐構(gòu)建一個健壯的計算函數(shù)了解了核心思路后我們需要將這些知識封裝成一個健壯、易用的MATLAB函數(shù)。這個函數(shù)不僅要能計算還要考慮實際應(yīng)用中的各種邊界情況和潛在錯誤。首先定義函數(shù)的目標(biāo)輸入一個包含至少8個元素的小時數(shù)據(jù)向量輸出其“最大8小時滑動平均值”。更完善一點我們還可以輸出整個8小時滑動平均序列以及最大值出現(xiàn)的位置結(jié)束小時。下面是一個綜合考慮了多種情況的函數(shù)實現(xiàn)示例function [max_8hr_avg, eight_hr_avg_series, max_idx] calc_max_8hr_moving_avg(hourly_concentration) % CALC_MAX_8HR_MOVING_AVG 計算小時濃度序列的最大8小時滑動平均值。 % % 輸入: % hourly_concentration - 數(shù)值向量按小時順序排列的濃度數(shù)據(jù)。 % % 輸出: % max_8hr_avg - 標(biāo)量最大8小時滑動平均值。 % eight_hr_avg_series - 向量與輸入等長的8小時滑動平均序列前7位為NaN。 % max_idx - 標(biāo)量最大值在 eight_hr_avg_series 中的索引結(jié)束小時。 % % 示例: % data randn(24,1)*10 50; % 模擬一天24小時數(shù)據(jù) % [maxVal, allAvg, idx] calc_max_8hr_moving_avg(data); % 1. 輸入驗證 if nargin 1 error(‘必須輸入小時濃度數(shù)據(jù)向量?！?; end if ~isvector(hourly_concentration) || ~isnumeric(hourly_concentration) error(‘輸入必須為數(shù)值向量?！?; end data hourly_concentration(:); % 強(qiáng)制轉(zhuǎn)換為列向量統(tǒng)一維度 N length(data); if N 8 error(‘輸入數(shù)據(jù)長度必須至少為8小時?！?; end % 2. 檢查數(shù)據(jù)有效性可選但很重要 % 假設(shè)無效數(shù)據(jù)如缺失值已被標(biāo)記為 NaN % 如果數(shù)據(jù)中包含NaNmovmean在默認(rèn)‘omitnan’模式下會忽略它們但這可能不符合某些標(biāo)準(zhǔn)。 % 這里我們采用嚴(yán)格模式窗口內(nèi)任何數(shù)據(jù)為NaN則結(jié)果也為NaN。 % 用戶可以根據(jù)需要修改此邏輯。 % 3. 核心計算使用 movmean window_size 8; % 關(guān)鍵參數(shù)[window_size-1, 0] 定義了非對稱窗口包含當(dāng)前點及前7個點。 % ‘Endpoints’, ‘fill’ 指定數(shù)據(jù)起始端不足窗口時用NaN填充。 % ‘omitnan’ 參數(shù)如果窗口內(nèi)存在NaN則計算結(jié)果為NaN。這是環(huán)境標(biāo)準(zhǔn)中常用的嚴(yán)格處理方式。 eight_hr_avg_series movmean(data, [window_size-1, 0], ‘Endpoints’, ‘fill’, ‘omitnan’); % 4. 尋找最大值忽略NaN [max_8hr_avg, max_idx] max(eight_hr_avg_series, ‘omitnan’); % 5. 處理全為NaN的特殊情況 if isnan(max_8hr_avg) max_8hr_avg NaN; max_idx NaN; warning(‘輸入的濃度數(shù)據(jù)可能全部為無效值NaN無法計算有效的滑動平均值。’); end end這個函數(shù)體現(xiàn)了幾個重要的工程化思考輸入驗證確保輸入是合法的數(shù)值向量且長度足夠。這是防止函數(shù)因意外輸入而崩潰的第一道防線。維度統(tǒng)一通過data(:)將輸入強(qiáng)制轉(zhuǎn)為列向量避免后續(xù)因行、列向量不同而導(dǎo)致的維度錯誤。NaN處理策略明確化在環(huán)境數(shù)據(jù)中缺失或無效的數(shù)據(jù)常以NaN表示。movmean的‘omitnan’選項會在計算窗口平均值時忽略NaN。但這需要特別注意如果一個8小時窗口內(nèi)只有4個有效數(shù)據(jù)movmean會計算這4個數(shù)據(jù)的平均值而不是8個。這不一定符合所有標(biāo)準(zhǔn)的規(guī)定。有些標(biāo)準(zhǔn)要求8小時內(nèi)必須有至少6個有效數(shù)據(jù)才計算平均值。因此在實際應(yīng)用中你可能需要先根據(jù)標(biāo)準(zhǔn)定義對原始數(shù)據(jù)中的NaN進(jìn)行預(yù)處理如插補或者編寫更復(fù)雜的邏輯來判斷窗口有效性。完整的輸出除了最大值還返回整個序列和索引便于用戶繪圖或進(jìn)一步分析最大值出現(xiàn)的時段。4. 實戰(zhàn)演練與深度避坑指南讓我們用一個更貼近現(xiàn)實的例子來演練。假設(shè)我們有一年8760小時的臭氧小時濃度模擬數(shù)據(jù)我們要計算每一天的“日最大8小時平均”。步驟1準(zhǔn)備測試數(shù)據(jù)% 生成模擬數(shù)據(jù)包含日周期、季節(jié)趨勢和隨機(jī)噪聲 hours_per_year 365*24; t (0:hours_per_year-1)‘; % 基礎(chǔ)水平 日周期白天高 年周期夏季高 噪聲 daily_cycle 30 * sin(2*pi*t/24 - pi/2) 30; % 峰值在下午 yearly_cycle 10 * sin(2*pi*t/(365.25*24)); noise randn(hours_per_year, 1) * 5; ozone_sim 20 daily_cycle yearly_cycle noise; ozone_sim max(ozone_sim, 0); % 濃度不為負(fù) % 故意插入一些缺失值NaN模擬真實數(shù)據(jù) missing_idx randi([1, hours_per_year], 100, 1); ozone_sim(missing_idx) NaN; % 將小時數(shù)據(jù)重塑為天數(shù) x 24小時的矩陣便于按天處理 ozone_matrix reshape(ozone_sim, 24, 365)‘; % 現(xiàn)在大小為 365 x 24步驟2逐日計算并可視化daily_max_8hr zeros(365, 1); daily_max_hour zeros(365, 1); % 記錄最大值結(jié)束的小時1-24 for day 1:365 hourly_data_day ozone_matrix(day, :); % 取出一行即一天24小時數(shù)據(jù) [max_val, ~, idx] calc_max_8hr_moving_avg(hourly_data_day’); daily_max_8hr(day) max_val; if ~isnan(idx) daily_max_hour(day) idx; else daily_max_hour(day) NaN; end end % 繪圖 figure(‘Position‘, [100, 100, 1200, 500]); subplot(2,1,1); plot(1:365, daily_max_8hr, ‘b-‘, ‘LineWidth‘, 1.5); xlabel(‘年積日‘); ylabel(‘日最大8小時平均臭氧濃度 (ppb)‘); title(‘模擬臭氧年變化日最大8小時平均值‘); grid on; subplot(2,1,2); scatter(1:365, daily_max_hour, 15, ‘filled‘); xlabel(‘年積日‘); ylabel(‘最大值出現(xiàn)的小時 (1-24)‘); title(‘最大值出現(xiàn)時刻分布‘); ylim([0 25]); grid on;步驟3關(guān)鍵避坑點與經(jīng)驗分享在實際操作中我踩過不少坑這里總結(jié)幾個最重要的時間序列的連續(xù)性與對齊這是最大的坑。我們的函數(shù)假設(shè)輸入數(shù)據(jù)是按小時嚴(yán)格連續(xù)排列的。但真實數(shù)據(jù)可能有時間戳。你必須確保你的hourly_concentration向量中的第i個元素確實對應(yīng)著第i個小時的濃度且中間沒有跳躍或缺失。如果數(shù)據(jù)有缺失小時直接計算會導(dǎo)致窗口錯位結(jié)果毫無意義。務(wù)必先進(jìn)行時間序列的重采樣或插值生成嚴(yán)格等間隔的連續(xù)序列。movmean的窗口定義movmean(data, [7, 0])和movmean(data, 8)天差地別。前者是非對稱窗口包含當(dāng)前點及前7點符合“向前追溯8小時”的定義。后者是對稱窗口包含當(dāng)前點、前3.5點和后3.5點MATLAB會自動處理為整數(shù)點這不符合環(huán)境標(biāo)準(zhǔn)。一定要根據(jù)你的業(yè)務(wù)需求精確選擇窗口模式。NaN的處理哲學(xué)如前所述‘omitnan’是雙刃劍。它讓你在數(shù)據(jù)有缺失時仍能得到一個數(shù)值但這個數(shù)值可能基于不完整的窗口。在撰寫報告或進(jìn)行達(dá)標(biāo)判斷時這可能導(dǎo)致錯誤結(jié)論。一個更穩(wěn)妥的做法是先定義一個“有效數(shù)據(jù)比例”閾值如6/875%。在計算每個窗口平均值前先判斷窗口內(nèi)非NaN數(shù)據(jù)的數(shù)量是否達(dá)標(biāo)若不達(dá)標(biāo)則直接給結(jié)果賦NaN。這需要自己用循環(huán)或movsum配合邏輯判斷來實現(xiàn)。計算“日最大”時的日期邊界問題環(huán)境標(biāo)準(zhǔn)中的“日最大8小時平均”通常是指“自然日”內(nèi)的最大值。但一個8小時窗口可能跨越兩天例如從第一天23點到第二天6點。嚴(yán)格來說這個跨日的8小時平均值應(yīng)該歸屬于第二天因為結(jié)束小時在第二天。在按天切片計算時如果你簡單地把每天0-23點單獨切片就會漏掉這些跨日窗口。正確的做法是在連續(xù)的長序列上先計算出所有8小時平均值然后根據(jù)每個平均值對應(yīng)的結(jié)束時間戳將其歸類到對應(yīng)的自然日中再在每個自然日里找最大值。這比簡單的按天循環(huán)更嚴(yán)謹(jǐn)。性能優(yōu)化對于超長序列如數(shù)十年每小時數(shù)據(jù)即使使用movmean一次性計算也可能內(nèi)存不足或速度較慢?？梢钥紤]使用tall array高數(shù)組或者將數(shù)據(jù)分塊處理。對于累積和法要注意數(shù)值精度問題當(dāng)數(shù)據(jù)量極大、數(shù)值跨度也大時cumsum可能導(dǎo)致浮點誤差累積但對于環(huán)境濃度數(shù)據(jù)通常問題不大。5. 進(jìn)階應(yīng)用擴(kuò)展到通用滑動窗口與性能對比我們的函數(shù)雖然解決了8小時的問題但我們可以很容易地將其泛化以計算任意窗口長度的滑動平均最大值。function [max_moving_avg, moving_avg_series, max_idx] calc_max_moving_avg(time_series, window_size, varargin) % CALC_MAX_MOVING_AVG 計算時間序列的指定窗口滑動平均最大值。 % 輸入: % time_series - 數(shù)值向量 % window_size - 正整數(shù)滑動窗口長度 % varargin - 可選參數(shù)對用于傳遞給 movmean如 ‘Endpoints’, ‘omitnan’ 等 % 輸出: (同上) % 輸入驗證 if nargin 2 error(‘必須輸入時間序列和窗口長度。’); end if ~isscalar(window_size) || window_size 1 || floor(window_size) ~ window_size error(‘窗口長度必須為正整數(shù)。’); end % ... (其他輸入驗證類似) % 核心計算允許自定義 movmean 參數(shù) if isempty(varargin) % 默認(rèn)參數(shù)非對稱窗口向前追溯起點填充NaN忽略NaN計算 moving_avg_series movmean(time_series, [window_size-1, 0], ‘Endpoints‘, ‘fill‘, ‘omitnan‘); else moving_avg_series movmean(time_series, [window_size-1, 0], varargin{:}); end % 尋找最大值 [max_moving_avg, max_idx] max(moving_avg_series, ‘omitnan‘); if isnan(max_moving_avg) max_moving_avg NaN; max_idx NaN; end end現(xiàn)在我們來對比一下之前提到的幾種實現(xiàn)方法在性能上的差異。我們用一個較大的數(shù)據(jù)集10萬個數(shù)據(jù)點來測試。% 性能測試 N 100000; test_data randn(N, 1); window_size 8; num_trials 100; % 運行次數(shù)取平均 % 方法1: movmean tic; for i 1:num_trials ma_mov movmean(test_data, [window_size-1, 0], ‘Endpoints‘, ‘fill‘); max_mov max(ma_mov, ‘omitnan‘); end time_movmean toc / num_trials; % 方法2: conv (valid模式) tic; for i 1:num_trials kernel ones(window_size, 1) / window_size; ma_conv conv(test_data, kernel, ‘valid‘); max_conv max(ma_conv); end time_conv toc / num_trials; % 方法3: 累積和法 tic; for i 1:num_trials cs [0; cumsum(test_data)]; ma_cumsum (cs(window_size1:end) - cs(1:end-window_size)) / window_size; max_cumsum max(ma_cumsum); end time_cumsum toc / num_trials; fprintf(‘性能對比 (窗口大小%d, 數(shù)據(jù)長度%d):\n‘, window_size, N); fprintf(‘ movmean: %.6f 秒\n‘, time_movmean); fprintf(‘ conv : %.6f 秒\n‘, time_conv); fprintf(‘ cumsum : %.6f 秒\n‘, time_cumsum);在我的測試環(huán)境MATLAB R2023b下結(jié)果通常是cumsum法最快conv法次之movmean稍慢但代碼最簡潔易讀。movmean作為內(nèi)置函數(shù)其優(yōu)勢在于功能豐富多種邊界處理、NaN處理選項并且代碼意圖一目了然。在大多數(shù)不是極端追求性能的場景下我強(qiáng)烈推薦使用movmean它的可讀性和維護(hù)性是最好的。只有當(dāng)處理海量數(shù)據(jù)且對速度有極致要求時才需要考慮手動實現(xiàn)累積和法。6. 在Simulink與實時系統(tǒng)中的應(yīng)用思考雖然本文主要討論在MATLAB命令窗口或腳本中的數(shù)據(jù)處理但滑動平均的概念在Simulink模型和實時嵌入式系統(tǒng)中也極其常見通常被稱為“移動平均濾波器”或“FIR均值濾波器”。在Simulink中你可以使用Moving Average模塊DSP System Toolbox或Discrete FIR Filter模塊將系數(shù)設(shè)為ones(1,8)/8來實現(xiàn)一個8點滑動平均濾波器。這對于處理來自傳感器的實時信號流非常有用。在實時C代碼實現(xiàn)中通常采用“循環(huán)緩沖區(qū)”來高效計算滑動平均避免每次都對窗口內(nèi)所有數(shù)據(jù)求和。其核心思路是維護(hù)一個窗口數(shù)據(jù)的和sum當(dāng)新數(shù)據(jù)x_new到來時減去從緩沖區(qū)中移出的最舊數(shù)據(jù)x_old再加上x_new然后計算平均值sum / N。這樣每次更新只需一次加法和一次減法復(fù)雜度是 O(1)非常適合單片機(jī)或嵌入式設(shè)備。// 偽代碼示例 float buffer[8]; // 循環(huán)緩沖區(qū) int index 0; // 當(dāng)前寫入位置 float sum 0; float update_moving_average(float new_sample) { // 減去即將被覆蓋的舊值 sum - buffer[index]; // 存入新值并加到總和里 buffer[index] new_sample; sum new_sample; // 更新索引 index (index 1) % 8; // 返回平均值 return sum / 8.0; }這種思路在MATLAB中也可以模擬對于處理流式數(shù)據(jù)很有啟發(fā)。但在處理完整的歷史數(shù)據(jù)集時我們更傾向于使用向量化的整體計算。最后無論是離線分析還是在線濾波理解滑動平均的頻率響應(yīng)特性都很重要。一個8點滑動平均濾波器是一個低通濾波器它會衰減高頻噪聲但同時也會使信號的快速變化變得“遲鈍”。其幅頻響應(yīng)是一個sinc函數(shù)在頻率為采樣頻率的1/8、2/8...等處會有零點。這意味著如果你的信號中有恰好是8小時倍數(shù)的周期成分會被完全濾除。在設(shè)計或解讀結(jié)果時需要意識到這個濾波器對數(shù)據(jù)本身特性的影響。回到我們最初的環(huán)境監(jiān)測例子計算“最大8小時平均”不僅僅是一個數(shù)學(xué)操作它背后是出于對人體健康風(fēng)險的考量——短期高暴露和長期平均暴露的影響是不同的。用MATLAB實現(xiàn)它是將業(yè)務(wù)規(guī)則轉(zhuǎn)化為可執(zhí)行代碼的典型過程。在這個過程中對細(xì)節(jié)的把握如窗口定義、NaN處理、時間對齊直接決定了結(jié)果的科學(xué)性和可靠性。希望這篇詳細(xì)的拆解能讓你下次面對類似“滑動平均”需求時不僅能寫出代碼更能理解每一行代碼背后的考量。