变分模态分解python

变分模态分解（VMD）是一种信号处理技术，可将信号分解成多个调制模态，每个模态都具有不同的频率和振幅。以下是使用Python进行VMD分解的基本步骤：

1.导入必要的库：numpy、scipy和matplotlib。

2.定义VMD函数：输入信号和参数（如模态数量、正则化参数等），并输出分解后的模态。

3.定义Hilbert-Huang变换（HHT）函数：将信号分解成固有模态和残差，以便进行VMD。

4.在主函数中，加载数据并调用VMD函数进行信号分解，并绘制每个模态及其相应的频谱。

下面是使用Python实现VMD的示例代码：

import numpy as np
from scipy.signal import hilbert
from matplotlib import pyplot as plt

def VMD(signal, alpha, tau, K, DC):
    # 参数设置
    N = len(signal)
    K = min(K, N)
    omega_k = np.zeros((K, N))
    u = signal
    u_hat = np.zeros((K, N))
    for k in range(K):
        # 做残差
        u_h = u - DC
        # 迭代分解
        for i in range(tau):
            u_hat[k, :] = hilbert(u_h)
            omega_k[k, :] = np.angle(u_hat[k, :])
            # 计算均值
            if k == K - 1:
                omega_mean = np.sum(omega_k, axis=0) / K
                omega_diff = omega_k - omega_mean
            else:
                omega_mean = np.sum(omega_k[:-1], axis=0) / (K - 1)
                omega_diff = omega_k[:-1] - omega_mean
            # 求出模态函数
            theta_k = np.exp(1j * omega_mean)
            u_h = u - DC - np.real(theta_k * u_hat[k, :])
        # 去讨论模态和直流分量
        DC += alpha * np.mean(u_h)
        u = u_h - alpha * np.mean(u_h)

    return u_hat

# 载入数据并进行VMD分解
signal = np.loadtxt('signal.txt')
K = 10  # 模态数量
alpha = 2500  # 正则化参数
tau = 10  # 迭代次数
DC = np.mean(signal)
u_hat = VMD(signal, alpha, tau, K, DC)

# 绘制每个模态及其频谱
plt.figure(figsize=(15, 10))
for i in range(K):
    # 绘制模态函数
    plt.subplot(K+1, 2, i*2+1)
    plt.plot(u_hat[i])
    plt.title(f'Mode {i+1}')
    plt.xlabel('Sample')
    plt.ylabel('Amplitude')
    # 绘制频谱
    plt.subplot(K+1, 2, i*2+2)
    freqs, psd = scipy.signal.welch(u_hat[i], fs=1000)
    plt.semilogy(freqs, psd)
    plt.title(f'PSD of Mode {i+1}')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power/Frequency (dB/Hz)')

# 绘制原始信号及其频谱
plt.subplot(K+1, 2, K*2+1)
plt.plot(signal)
plt.title('Original Signal')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.subplot(K+1, 2, K*2+2)
freqs, psd = scipy.signal.welch(signal, fs=1000)
plt.semilogy(freqs, psd)
plt.title('PSD of Original Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power/Frequency (dB/Hz)')
plt.tight_layout()
plt.show()