matlab程序实现图片pod本征正交分解的模态能量分布,累积能量分布,时间系数利萨如图,时间系数时序图和频谱图
img_dir = 'snapshots/'; file_list = dir(fullfile(img_dir,'*.png')); snapshots = []; for k = 1:length(file_list) img = imresize(imread(fullfile(img_dir,file_list(k).name)), [256 256]); snapshots(:,k) = double(img(:)); % 强制拉成列向量堆叠 end这里有个坑:不同尺寸图片得提前统一大小,imresize直接暴力缩放。数据矩阵每列代表一个时间点的全场信息,行数就是像素总数。
核心POD分解其实就是玩SVD:
[U,S,V] = svd(snapshots, 'econ'); lambda = diag(S).^2; % 特征值平方才是能量 energy_ratio = lambda / sum(lambda); cum_energy = cumsum(energy_ratio);U矩阵的列就是POD模态,S对角阵存奇异值。这里注意svd默认全分解,实测用'econ'选项能省内存。能量占比计算别手滑忘了归一化,cumsum直接给累积曲线。
matlab程序实现图片pod本征正交分解的模态能量分布,累积能量分布,时间系数利萨如图,时间系数时序图和频谱图
画能量分布建议用semilogy,一眼看出模态重要性:
figure('Position',[200 200 800 300]) subplot(1,2,1) plot(energy_ratio(1:20),'bo-','LineWidth',1.5) xlabel('Mode Number'); title('Energy Distribution') subplot(1,2,2) plot(cum_energy(1:20),'rs--','LineWidth',1.5) hold on; yline(0.95,'--'); xlabel('Mode Number'); title('Cumulative Energy')通常前5%的模态扛了95%以上的能量,这就是POD牛逼之处——用少量模态近似全场。
时间系数这块,V矩阵乘以S就是时间演化:
time_coeffs = S*V'; t = 1:size(time_coeffs,2); % 假设时间步长均匀 figure plot(time_coeffs(1,:), time_coeffs(2,:), '.') % 利萨如图 xlabel('Mode 1'); ylabel('Mode 2')如果看到闭合曲线,说明存在周期性运动。再补个时序图:
figure subplot(2,1,1) plot(t, time_coeffs(1,:), 'b') xlabel('Time'); ylabel('Amplitude') subplot(2,1,2) Fs = 1000; % 采样频率自己按实际情况改 L = length(t); Y = fft(time_coeffs(1,:)); P2 = abs(Y/L); P1 = P2(1:L/2+1); f = Fs*(0:(L/2))/L; plot(f, P1) xlabel('f (Hz)'); ylabel('|Amplitude|')FFT前注意去均值,否则零频分量巨高。频谱图里冒尖儿的频率对应主导周期,这对流场分析特别有用。
最后说个经验:实际计算时快照别超过2000个,否则SVD算到地老天荒。内存不够的话改用快照POD法,转置矩阵再分解,自己查文献吧。