利用python 实现dmd,且进行特征圆的提取，频谱图的提取

下面是一个完整的DMD实现示例，包括特征圆和频谱图的提取：

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import svd

def DMD(X, r):
    # X: 输入数据，为一个(m, n)的矩阵，m为数据点数，n为时间点数
    # r: 分解时取前r个奇异值
    
    # 构造数据矩阵X1和X2
    X1 = X[:, :-1]
    X2 = X[:, 1:]
    
    # 奇异值分解
    U, S, V = svd(X1, full_matrices=False)
    
    # 取前r个奇异值，构造S+矩阵
    Sr = np.diag(S[:r])
    Up = U[:, :r]
    Vp = V[:r, :]
    
    # 计算A矩阵
    A = Up.T @ X2 @ Vp @ np.linalg.inv(Sr)
    
    # 计算特征值和特征向量
    eigvals, eigvecs = np.linalg.eig(A)
    
    # 计算DMD模式
    modes = X2 @ Vp @ np.linalg.inv(Sr) @ eigvecs
    
    return modes, eigvals

# 生成数据
t = np.linspace(0, 2*np.pi, 200)
x1 = np.exp(1j*t)
x2 = np.exp(2j*t)
X = np.vstack([x1, x2])

# DMD分解
r = 2
modes, eigvals = DMD(X, r)

# 提取特征圆
center = np.log(eigvals) / (2*np.pi)
radius = np.abs(np.log(eigvals)) / (2*np.pi)

# 绘制特征圆
fig, ax = plt.subplots()
ax.scatter(center.real, center.imag, s=100, c='r', marker='x', label='feature circle')

# 提取频谱图
freqs = np.angle(eigvals) / (2*np.pi)
amplitudes = np.abs(np.diag(eigvals))
spectrum = amplitudes.reshape(-1, 1) @ np.exp(2j*np.pi*freqs).reshape(1, -1)

# 绘制频谱图
fig, ax = plt.subplots()
ax.imshow(np.log10(1 + spectrum), cmap='jet', aspect='auto', origin='lower', extent=[0, 1, 0, 1])
ax.set_xlabel('Frequency')
ax.set_ylabel('Time')

plt.show()

这个DMD实现生成两个信号，然后进行DMD分解，并提取特征圆和频谱图。特征圆可以用于判断DMD分解结果的稳定性和特征模式的周期性，频谱图可以用于分析信号的频域特征。