matlab动态模态分解程序

%这是一个动态模态分解matlab程序，X和Y已知，对此数据进行分析，我想将随时间衰减很快，且振幅（能量）大的模态都减掉，应该怎么写程序？
（实在不好意思修改了一下题目）

[Dd,b,Phi,Time_DMD,Energy] = DMD_CLASS(X,Y);
function [Dd,b,Phi,Time_DMD,Energy] = DMD_CLASS(X,Y)
%DMD分解函数
%输入：
%X，Y，DMD分解的数据矩阵。X是data(1,1:(end-1))，Y是data(:,2:end)。
%
%输出：
%Dd，特征根
%b，初始状态
%Phi，DMD分解的模态
%Time_DMD，DMD还原信号所用到的时间项
%Energy，每个模态对应的能量，从大到小排序

N=size(X,2);
%1 SVD分解
[U,S,V] = svd(X,'econ');%奇异值分解的精简分解
%删除奇异值约等于0的模态，防止计算发散
Sd = diag(S);
r = sum(Sd>1e-4);%可以将条件改成1e-6，1e-5等
U = U(:,1:r);
S = S(1:r,1:r);
V = V(:,1:r);
%2 求解矩阵A
A = U'*(Y*V/S);
%3 求矩阵A的特征值和特征向量
[Om,D] = eig(A);
Dd = diag(D);%求特征值转为向量形式
%4 求DMD模态
Phi = Y*V/S*Om;
%5 求解初始状态b
b = Phi\X(:,1);
%6 对模态进行排序
Q = Dd.^(0:N-1);%计算出范德蒙矩阵
Time_DMD = b.*Q;%旧版本matlab可以用右边形式替换Time_DMD=(b*ones(1,N)).*Q;

% %对信号初始振幅排序，速度快，输出能量可以近似用振幅平方代替
%[b,Ie] = sort(abs(b),'descend');
%Energy = abs(b).*abs(b);

%求出所有模态的信号，并计算能量。这是按照能量排序，和上面3行不一样。毕竟DMD模态的排序方法不唯一。
Energy = zeros( size(Phi,2) , 1 );
for k = 1 : size(Phi,2)
Uxt_DMD_k = real(Phi(:,k) * Time_DMD(k,:));
E_k = sum( sum( Uxt_DMD_k.^2 ) );
Energy(k) = E_k;
end
[Energy,Ie] = sort(Energy,'descend');%对每个模态的能量进行排序

%按顺序输出特征值、初始状态b、模态Phi
Dd = Dd(Ie);
b = b(Ie);
Phi = Phi(:,Ie);
Time_DMD = Time_DMD(Ie,:);
end

这是一个组合函数吧 , 你少给我东西了吧?

要减去随时间衰减的模式，需要识别并删除衰减率为负或较小的模式。可以通过添加负或小衰减率检查来修改DMD_CLASS函数，并将这些模式设置为零。示例：

function [Dd,b,Phi,Time_DMD,Energy] = DMD_CLASS(X,Y,delta)
%DMD分解函数
%输入：%X，Y，DMD分解的数据矩阵。X是data(1,1:(end-1))，Y是data(:,2:end)。%
%delta, 阈值，用于判断需要保留的模态的衰减率，即如果模态的特征值小于delta，则该模态被认为衰减严重，需要被移除。
%%输出：%Dd，特征根
%b，初始状态
%Phi，DMD分解的模态
%Time_DMD，DMD还原信号所用到的时间项
%Energy，每个模态对应的能量，从大到小排序

N=size(X,2);

%1 SVD分解
[U,S,V] = svd(X,'econ');%奇异值分解的精简分解

%删除奇异值约等于0的模态，防止计算发散
Sd = diag(S);
r = sum(Sd>1e-4);%可以将条件改成1e-6，1e-5等
U = U(:,1:r);
S = S(1:r,1:r);
V = V(:,1:r);

%2 求解矩阵A
A = U'*(Y*V/S);

%3 求矩阵A的特征值和特征向量
[Om,D] = eig(A);
Dd = diag(D);%求特征值转为向量形式

%4 求DMD模态
Phi = Y*V/S*Om;

%5 求解初始状态
bb = Phi\X(:,1);

%6 对模态进行排序
Q = Dd.^(0:N-1);%计算出范德蒙矩阵
Time_DMD = b.*Q;%旧版本matlab可以用右边形式替换Time_DMD=(b*ones(1,N)).*Q;

%对信号初始振幅排序，速度快，输出能量可以近似用振幅平方代替
[b,Ie] = sort(abs(b),'descend');
Energy = abs(b).*abs(b);%求出所有模态的信号，并计算能量。这是按照能量排序，和上面3行不一样。毕竟DMD模态的排序方法不唯一。
% Energy = zeros( size(Phi,2) , 1 ); % use this if you want to compute energy using reconstructed signals
for k = 1 : size(Phi,2)
    % skip modes with negative or small decay rates
    if abs(Dd(k)) < delta
        Phi(:,k) = 0;
        continue;
    end
    
    Uxt_DMD_k = real(Phi(:,k) * Time_DMD(k,:));
    % E_k = sum( sum( Uxt_DMD_k.^2 ) ); % use this if you want to compute energy using reconstructed signals
    Energy(k) = abs(Dd(k)); % set energy to decay rate
end

[Energy,Ie] = sort(Energy,'descend');%对每个模态的能量进行排序

%按顺序输出特征值、初始状态b、模态Phi
Dd = Dd(Ie);
b = b(Ie);
Phi = Phi(:,Ie);
Time_DMD = Time_DMD(Ie,:);
end

该修改包括一个附加的输入参数“delta”，它设置了衰减率的阈值，低于该阈值应去除模式。在计算每个模式能量的循环中，对衰减率（“abs（Dd（k））＜delta”）进行检查，如果它小于阈值，则将模式设置为零。

要使用这个修改后的函数，只需用数据矩阵“X”、“Y”和所需的衰减率阈值“delta”来调用它。输出将与原始的“DMD_CLASS”函数相同，只是衰减过快的模式已被删除。

可以这样做，借鉴一下这个：https://zhuanlan.zhihu.com/p/635043276

参考代码


function [Dd,b,Phi,Time_DMD,Energy] = DMD_CLASS(X,Y, energy_thresh)

N = size(X,2);

[U,S,V] = svd(X,'econ');

Sd = diag(S);
r = sum(Sd>1e-4);
U = U(:,1:r);
S = S(1:r,1:r);
V = V(:,1:r);

A = U'*(Y*V/S);

[Om,D] = eig(A);
Dd = diag(D);

Phi = Y*V/S*Om;

b = Phi\X(:,1);

Q = Dd.^(0:N-1);
Time_DMD = b.*Q;

Energy = zeros(size(Phi,2), 1);
for k = 1:size(Phi,2)
    Uxt_DMD_k = real(Phi(:,k) * Time_DMD(k,:));
    %计算模态在最后一个时间步长的振幅大小
    b(k) = Uxt_DMD_k(end);
    %计算模态的能量
    E_k = sum(sum(Uxt_DMD_k.^2));
    Energy(k) = E_k;
end

%找到能量小于阈值的模态并删除
energy_pct = (Energy ./ sum(Energy)) * 100;
mode_inds_to_keep = find(energy_pct >= energy_thresh);
Dd = Dd(mode_inds_to_keep);
Phi = Phi(:,mode_inds_to_keep);
b = b(mode_inds_to_keep);
Time_DMD = Time_DMD(mode_inds_to_keep,:);

%按照能量大小对模态排序
[Energy, sort_inds] = sort(Energy, 'descend');
Dd = Dd(sort_inds);
Phi = Phi(:,sort_inds);
b = b(sort_inds);
Time_DMD = Time_DMD(sort_inds,:);

end

可以参考下

%DMD实验2
clear
close all

x=0:0.01:5;m=length(x);
t=0:0.1:5;N=length(t);
Fs=1/(t(2)-t(1));
[X,T]=meshgrid(x,t);

%信号
U_tx=1.2*exp(-0.5*T).*sin(2*pi*(X+0*T))+0.8*exp(0.3*T).*sin(2*pi*(3*X+4*T));
%p的格式为时间*空间的矩阵
U_xt=U_tx';
Up=pinv(U_xt(:,1:50));%求解X的伪逆
%求解矩阵A
A=U_xt(:,2:51)*Up;

[V,D] = eig(A);
D=diag(D);
[~,ind] = sort(abs(D),'descend');
D=D(ind);V=V(:,ind);
D_abs=abs(D);
Dd=D./D_abs;
%频率
f=Fs*angle(Dd)./(2*pi);

%图1，绘制特征根图谱
figure(1)
scatter(real(Dd(100:-1:1)),imag(Dd(100:-1:1)),30,-abs(D(100:-1:1)),'filled')
axis equal
box on

%图2，绘制频率和衰减图
figure(2)
wa=log(D)*Fs;%利用特征根计算衰减率
scatter(real(wa),imag(wa)/2/pi,30,-D_abs,'filled')
xlabel('衰减率σ');ylabel('频率w')
ylim([-6,6]);xlim([-1,1])
hold on
plot([0,0],ylim,'b--')
plot(xlim,[0,0],'b--')
box on

D模态对应的能量较小，随时间衰减快，可以通过计算每个DMD模态对应的能量，然后去掉能量较高的模态来减少振幅较高的数据。代码中的变量Energy就是每个模态对应的能量，可以通过对能量从大到小排序，确定需要去掉的模态数目，将这些模态对应的数据从原始数据中剔除即可。

具体实现代码如下：

Energy = sort(abs((diag(S)V')b).^2/(prod(abs(Sd))trace(YY')),'descend'); %计算所有模态对应的能量,从大到小排序 total_energy = sum(Energy); %计算所有模态对应的总能量 energy_ratio = cumsum(Energy)/total_energy; %计算累计能量贡献率 mode_num = find(energy_ratio>=0.95,1); %确定需要去掉的模态数目 modes_to_remove = 1:mode_num; %需要去掉的模态 keep_modes = mode_num+1:length(Energy); %需要保留的模态 DMD_modes(:,modes_to_remove) = []; %剔除需要去掉的模态 DMD_modes_energy = Energy(keep_modes); %保留的模态对应的能量

根据您提供的MATLAB代码，您想要对给定的数据进行动态模态分解，并减去随时间衰减的模态。以下是编写程序的步骤：

调用函数：调用DMD_CLASS函数，并传入数据矩阵X和Y。
SVD分解：对数据矩阵X进行奇异值分解（SVD），得到矩阵U、S和V。根据奇异值的大小，选择一个阈值（例如1e-4），确定需要保留的奇异值个数r，并截取U、S和V的前r列。
计算矩阵A：根据公式 A = U'_(Y_V/S)，计算矩阵A。
求解特征值和特征向量：对矩阵A进行特征值分解，得到特征值D和特征向量Om。
计算DMD模态：根据公式 Phi = Y_V/S_Om，计算DMD模态矩阵Phi。
求解初始状态b：根据公式 b = Phi\X(:,1)，求解初始状态b。
对模态进行排序：计算范德蒙矩阵Q和Time_DMD，然后根据能量排序。
输出结果：输出特征根Dd、初始状态b、DMD模态Phi、DMD还原信号所用的时间项Time_DMD以及每个模态对应的能量Energy。

根据您的需求，如果您想将随时间衰减的模态减掉，可以在计算DMD还原信号时，根据需要减去衰减较快的模态。

注意：以上步骤是根据您提供的代码进行理解和解释的，确保数据矩阵X和Y的格式正确，并适用于您的具体问题。

根据你的描述，你希望将随时间衰减很快且振幅（能量）大的模态减掉。你可以在函数 DMD_CLASS 中的最后几行进行修改，按照以下步骤编写程序：

将以下两行代码注释掉，因为你不再需要按初始振幅排序和计算能量：


% [b,Ie] = sort(abs(b),'descend');
% Energy = abs(b).*abs(b);

将以下代码取消注释，以按照能量排序模态：

[Energy,Ie] = sort(Energy,'descend'); % 对每个模态的能量进行排序

修改模态的输出，只保留衰减较慢且能量较小的模态：

Dd = Dd(Ie);
b = b(Ie);
Phi = Phi(:,Ie);
Time_DMD = Time_DMD(Ie,:);

通过这些修改，模态将按照能量从大到小排序，并且你可以选择保留衰减较慢且能量较小的模态。你可以根据需要修改代码，例如只保留能量占总能量一定比例的模态。

对于动态模态分解（DMD）程序，可以通过对DMD模态的能量进行排序，然后选择能量较大的模态进行保留，而将能量较小的模态进行截断。这样可以达到减掉随时间衰减很快、且振幅（能量）小的模态的效果。

具体实现步骤如下：

运行DMD程序，得到DMD模态和对应的能量：
[Phi,omega,lambda,b,X_dmd,S] = DMD(X,Y,r);
其中， Phi 是DMD模态， omega 是DMD模态的频率， lambda 是DMD模态的特征值， b 是DMD重构系数， X_dmd 是DMD重构结果， S 是DMD模态的能量。
对能量进行排序：
[sorted_S, index] = sort(S, 'descend');
这里使用 sort 函数将能量从大到小进行排序，并返回排序后的能量和对应的索引。
选择能量较大的模态进行保留：
k = 10; % 选择前10个能量较大的模态
Phi_k = Phi(:, index(1:k));
omega_k = omega(index(1:k));
lambda_k = lambda(index(1:k));
b_k = b(index(1:k));
这里选择前10个能量较大的模态进行保留，可以根据实际情况进行调整。
对DMD进行重构：
X_dmd_k = Phi_k * diag(b_k) * exp(omega_k * t);
这里只使用了保留的模态进行重构，得到新的重构结果 X_dmd_k 。

完整的代码如下：

% 运行DMD程序，得到DMD模态和对应的能量
[Phi,omega,lambda,b,X_dmd,S] = DMD(X,Y,r);
 % 对能量进行排序
[sorted_S, index] = sort(S, 'descend');
 % 选择前10个能量较大的模态进行保留
k = 10;
Phi_k = Phi(:, index(1:k));
omega_k = omega(index(1:k));
lambda_k = lambda(index(1:k));
b_k = b(index(1:k));
 % 对DMD进行重构
X_dmd_k = Phi_k * diag(b_k) * exp(omega_k * t);