该算法是一个综合性的信号处理与分析框架，专门设计用于提取和分析复杂信号中的多尺度特征与多重分形特性。算法的核心思想是通过变分模态分解VMD将原始复杂信号自适应分解为多个相对简单的本征模态函数IMF，然后对每个模态进行多重分形分析MFA，以揭示信号在不同频率尺度上的复杂结构和特性。具体而言，算法首先生成或加载一个具有多重分形特性的信号，如多重分形随机游走MRW，然后使用VMD方法将信号分解为多个频率成分不同的模态，每个模态代表信号在不同频率范围内的振荡成分。接着，算法对每个模态进行小波变换和p-leaders计算，进而执行多重分形分析，包括计算标度函数ζ(q)、累积量和多重分形谱D(h)。通过这些分析，算法能够量化每个模态的多重分形特性，揭示信号在不同尺度上的复杂行为。最后，算法通过多种可视化手段展示分解结果和分形特性，包括时域波形、频谱分析、中心频率演化、功率谱密度和模态间相关性等，为用户提供全面的信号特征分析结果。

开始
│
├─ 数据准备
│   ├─ 生成多重分形随机游走(MRW)信号
│   └─ 可视化原始信号
│
├─ 完整信号多重分形分析
│   ├─ 小波变换(db3小波)
│   ├─ 计算p-leaders (p=2)
│   ├─ 在指定尺度范围内积分
│   ├─ 多重分形分析(计算ζ(q)、累积量、D(h))
│   └─ 可视化分析结果
│
├─ 变分模态分解(VMD)
│   ├─ 设置VMD参数(alpha, tau, K, DC, init, tol)
│   ├─ 执行VMD分解
│   ├─ 按频率排序模态
│   └─ 可视化分解结果
│
├─ 模态特性分析
│   ├─ 绘制所有模态时域波形
│   ├─ 分别绘制每个模态
│   ├─ 绘制频谱分解结果(对数坐标)
│   ├─ 绘制中心频率演化过程
│   ├─ 计算功率谱密度(Welch方法)
│   └─ 计算模态间相关系数
│
└─ 模态多重分形分析
    ├─ 对每个模态进行小波变换
    ├─ 计算p-leaders
    ├─ 在指定尺度范围内积分
    ├─ 多重分形分析
    └─ 可视化每个模态的分形特性

算法步骤详解

1. 数据准备阶段：生成多重分形随机游走(MRW)信号作为示例数据，MRW是一种常用的多重分形过程模拟方法，能够产生具有复杂尺度特性的时间序列数据，可视化原始信号以直观了解其基本特征。
2. 完整信号多重分形分析阶段：对完整信号进行多重分形分析，使用db3小波进行小波变换，计算p-leaders(p=2)以捕获信号的局部正则性，在指定尺度范围内进行积分以聚焦特定尺度范围，执行多重分形分析计算标度函数ζ(q)、累积量和多重分形谱D(h)，可视化分析结果以展示信号的全局多重分形特性。
3. 变分模态分解阶段：设置VMD分解参数包括带宽约束、噪声容忍度、模态数量等，执行VMD分解将原始信号分解为多个本征模态函数，按最终频率排序模态以便分析，可视化分解结果包括所有模态的时域波形和每个模态的单独显示。
4. 模态特性分析阶段：绘制所有模态的时域波形以直观比较各模态的时域特征，分别绘制每个模态以详细分析单个模态的特性，绘制频谱分解结果在对数坐标下展示各模态的频域分布，绘制中心频率演化过程展示VMD迭代过程中各模态频率的变化，计算功率谱密度使用Welch方法量化各模态的频域能量分布，计算模态间相关系数评估各模态在时域上的独立性。
5. 模态多重分形分析阶段：对每个VMD模态分别进行小波变换以提取尺度相关信息，计算p-leaders以量化每个模态的局部正则性，在指定尺度范围内进行积分以聚焦分析范围，执行多重分形分析计算每个模态的标度函数、累积量和多重分形谱，可视化每个模态的分形特性以比较不同模态的多重分形行为差异。

# 导入必要的库
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from pymultifracs.utils import build_q_log
from pymultifracs.simul import mrw
from pymultifracs.wavelet import wavelet_analysis
from pymultifracs.mf_analysis import mfa
from sktime.libs.vmdpy import VMD
from scipy.signal import welch

# 设置绘图风格
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("husl")

# 生成多重分形随机游走(MRW)信号作为示例数据
# MRW是一种常用的多重分形过程模拟方法
nb_generated_series = 1  # 生成1个时间序列
X = mrw(shape=(2**15, nb_generated_series), H=0.8, lam=np.sqrt(.05), L=2**15)

# 绘制生成的MRW信号
plt.figure(figsize=(12, 5))
plt.plot(X)
plt.title('Multifractal Random Walk (H=0.8, λ=√0.05)')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

# 定义多重分形分析的尺度范围
scaling_ranges = [[6, 11]]

# 对完整信号进行多重分形分析
# 使用db3小波进行小波分析
WT = wavelet_analysis(X, wt_name='db3')
# 计算p-leaders (p=2)
WTpL = WT.get_leaders(p_exp=2)
# 在指定尺度范围内进行积分
WTpL = WTpL.auto_integrate(scaling_ranges)
# 执行多重分形分析，使用对数间隔的q值
pwt = mfa(WTpL, scaling_ranges, weighted='Nj', q=build_q_log(.1, 5, 20))

# 绘制完整信号的标度函数ζ(q)
plt.figure(figsize=(8, 5))
pwt.structure.plot_scaling()
plt.title("ζ(q) for The Original Signal")
plt.xlabel("q")
plt.ylabel("ζ(q)")
plt.tight_layout()
plt.show()

# 绘制完整信号的累积量
plt.figure(figsize=(8, 5))
pwt.cumulants.plot()
plt.title("Cumulants for The Original Signal")
plt.xlabel("Scale j")
plt.ylabel("Cumulant Value")
plt.tight_layout()
plt.show()

# 绘制完整信号的多重分形谱D(h)
plt.figure(figsize=(8, 5))
pwt.spectrum.plot()
plt.title("D(h) Spectrum for the Original Signal")
plt.xlabel("Singularity Exponent h")
plt.ylabel("Dimension D(h)")
plt.tight_layout()
plt.show()

# 提取信号进行VMD分解
signal = X[:, 0]  # 提取一维信号
N = len(signal)   # 信号长度

# 计算信号的傅里叶变换
frequencies = np.fft.fftshift(np.fft.fftfreq(N, d=1/N)) * 2 * np.pi
signal_fft = np.fft.fftshift(np.fft.fft(signal))

# 设置VMD参数
alpha = 2000  # 带宽约束参数
tau = 0.      # 噪声容忍度
K = 3         # 分解模态数
DC = 0        # 不包含直流分量
init = 1      # 初始化方式(均匀分布)
tol = 1e-7    # 收敛容差

# 执行VMD分解
u, u_hat, omega = VMD(signal, alpha, tau, K, DC, init, tol)

# 按最终频率排序模态
sortIndex = np.argsort(omega[-1, :])
omega = omega[:, sortIndex]
u_hat = u_hat[:, sortIndex]
u = u[sortIndex, :]

# 绘制所有模态的时域波形
plt.figure(figsize=(12, 8))
for k in range(K):
    plt.plot(u[k, :], label=f"Mode {k+1}")
plt.title("Decomposed Modes from MRW Signal")
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 分别绘制每个模态
plt.figure(figsize=(12, 10))
for k in range(K):
    plt.subplot(K, 1, k + 1)
    plt.plot(u[k, :])
    plt.title(f"Decomposed Mode {k+1}")
    plt.xlabel("Time")
    plt.ylabel("Amplitude")
    plt.xlim(0, N)
plt.tight_layout()
plt.show()

# 绘制频谱分解结果(对数坐标)
plt.figure(figsize=(10, 6))
plt.loglog(frequencies[N//2:], np.abs(signal_fft[N//2:]), 'k:', label='Original')
for k in range(K):
    plt.loglog(frequencies[N//2:], np.abs(u_hat[N//2:, k]), label=f'Mode {k+1}')
plt.xlim([frequencies[N//2], frequencies[-1]])
plt.xlabel("Frequency (rad/s)")
plt.ylabel("Amplitude (log scale)")
plt.title("Spectral Decomposition")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 绘制中心频率演化过程
fs = 1 / N  # 采样频率近似
plt.figure(figsize=(10, 6))
for k in range(K):
    plt.semilogx(2 * np.pi / fs * omega[:, k], 
                 np.arange(1, omega.shape[0] + 1),
                 label=f"Mode {k+1}")
plt.title("Evolution of Center Frequencies ωₖ (VMD iterations)")
plt.xlabel("Frequency (rad/s)")
plt.ylabel("Iteration")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend()
plt.tight_layout()
plt.show()

# 使用Welch方法计算功率谱密度
plt.figure(figsize=(10, 6))
for k in range(K):
    f_welch, Pxx = welch(u[k, :], fs=1.0, nperseg=1024)
    plt.semilogy(f_welch, Pxx, label=f"Mode {k+1}")
plt.title("Mode Spectra via Welch PSD")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power / Hz")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 计算模态间的相关系数
df = pd.DataFrame(u.T, columns=[f'Mode {i+1}' for i in range(K)])
corr_matrix = df.corr()

plt.figure(figsize=(8, 6))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0)
plt.title("Correlation Between Modes (Time Domain)")
plt.tight_layout()
plt.show()

# 对每个VMD模态进行多重分形分析
for k in range(K):
    # 小波分析
    WT = wavelet_analysis(u[k, :], wt_name='db3')
    
    # 计算p-leaders
    WTpL = WT.get_leaders(p_exp=2)
    
    # 在指定尺度范围内积分
    WTpL = WTpL.auto_integrate(scaling_ranges)
    
    # 多重分形分析
    pwt = mfa(WTpL, scaling_ranges, weighted='Nj', q=build_q_log(.1, 5, 20))

    # 绘制标度函数ζ(q)
    plt.figure(figsize=(8, 5))
    pwt.structure.plot_scaling()
    plt.title(f"ζ(q) for Mode {k+1}")
    plt.xlabel("q")
    plt.ylabel("ζ(q)")
    plt.tight_layout()
    plt.show()

    # 绘制累积量
    plt.figure(figsize=(8, 5))
    pwt.cumulants.plot()
    plt.title(f"Cumulants for Mode {k+1}")
    plt.xlabel("Scale j")
    plt.ylabel("Cumulant Value")
    plt.tight_layout()
    plt.show()

    # 绘制多重分形谱D(h)
    plt.figure(figsize=(8, 5))
    pwt.spectrum.plot()
    plt.title(f"D(h) Spectrum for Mode {k+1}")
    plt.xlabel("Singularity Exponent h")
    plt.ylabel("Dimension D(h)")
    plt.tight_layout()
    plt.show()

知乎学术咨询（哥廷根数学学派）

工学博士，担任《Mechanical System and Signal Processing》审稿专家，担任

《中国电机工程学报》《控制与决策》，《系统工程与电子技术》，《电力系统保护与控制》，《宇航学报》等EI期刊审稿专家，担任《计算机科学》，《电子器件》 ，《现代制造过程》 ，《电源学报》，《船舶工程》 ，《轴承》 ，《工矿自动化》 ，《重庆理工大学学报》 ，《噪声与振动控制》 ，《机械传动》 ，《机械强度》 ，《机械科学与技术》 ，《机床与液压》，《声学技术》，《应用声学》等中文核心审稿专家。

擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。