开发工具:
文件大小: 7kb
下载次数: 0
上传时间: 2019-06-02
详细说明： import numpy as np from numba import jit jit def VMD(signal, fs, T, alpha, tau, K, DC, init, tol): f_mirror = np.zeros((2*T)) for i in range(T//2): f_mirror[i] = signal[T//2-1-i] for k in range(T//2, 3*T//2): f_mirror[k] = signal[k - T // 2] for l in range(3*T//2, 2*T): f_mirror[l] = signal[5 *T//2-l-1] f = f_mirror T = len(f) t = np.arange(1, T+1, 1)/T freqs = t - 0.5 - 1.0 / T freqs = freqs + 0.j N = 500 Alpha = alpha * np.ones((1, K)) Alpha = Alpha + 0.j f_hat = np.fft.fftshift(np.fft.fft(f)) f_hat_plus = f_hat for l in range(0, T//2): f_hat_plus[l] = 0 u_hat_plus = np.zeros((N, T, K)) u_hat_plus = u_hat_plus + 0.j omega_plus = np.zeros((N, K)) # 生成N行K列的数组，元素为0， 500*3 omega_plus = omega_plus + 0.j if init == 1: for i in range(0, K): omega_plus[0][i] = 0.5/K*i elif init == 2: omega_plus[0, :] = np.sort(np.exp(np.log(fs)) + (np.log(0.5)-np.log(fs))*np.random.rand(0, K)) else: omega_plus[0, :] = 0 if DC: omega_plus[0, 0] = 0 else: pass # 从空双变量开始 lambda_hat = np.zeros((N, T)) # lambda_hat = lambda_hat + 0.j # 其他参数初始化 uDiff = tol + 2.77555756e-17 # 更新步长 n = 1 # 循环计数 sum_uk = 0.j # 器材累加 while uDiff > tol and n < N : # 表示没有收敛或低于迭代次数，其中,n和K都从1开始 # 更新第一阶模态累加器 k = 0 sum_uk = u_hat_plus[n-1, :, K-1] + sum_uk - u_hat_plus[n-1, :, 0] # sum_uk此时还是一维数组 # print(sum_uk) # ------------------------ ----检验专用----------------------- # 通过残差维纳滤波器更新第一模态的频谱 u_hat_plus[n, :, k] = (f_hat_plus - sum_uk - lambda_hat[n - 1, :] / 2.0) / (1 + Alpha[0, k]*(freqs - omega_plus[n-1, k])**2) if ~DC: omega_plus[n, k] = np.dot(freqs[T//2: T], ((abs(u_hat_plus[n, T//2:T, k])).reshape(T//2, 1))**2) / sum((abs(u_hat_plus[n, T//2:T, k]))**2) for k in range(1, K): # 累加 sum_uk = u_hat_plus[n, :, k-1] + sum_uk - u_hat_plus[n-1, :, k] # print(sum_uk[1579], sum_uk[1579]) # -----------------------------检验专用----------------------- # 模态频谱 u_hat_plus[n, :, k] = (f_hat_plus - sum_uk - lambda_hat[n-1, :] / 2.0)/(1 + Alpha[0, k]*(freqs - omega_plus[n-1, k])**2) # print(u_hat_plus[n, 1001, k], u_hat_plus[n, 1001, k]) # 中心频率 omega_plus[n, k] = np.dot(freqs[T//2: T], (abs(u_hat_plus[n, T//2:T, k])**2).reshape(T//2, 1))/sum(abs(u_hat_plus[n, T//2:T, k])**2) # 双重上升 lambda_hat[n, :] = lambda_hat[n-1, :] + tau*(np.sum(u_hat_plus[n, :, :], 1) - f_hat_plus) n = n + 1 uDiff = 2.77555756e-17 for i in range(0, K): uDiff = uDiff + 1/T * np.dot((u_hat_plus[n-1, :, i] - u_hat_plus[n - 2, :, i]), (np.conj((u_hat_plus[n-1, :, i] - u_hat_plus[n-2, :, i]))).reshape(T, 1)) uDiff = abs(uDiff) N = min(N, n) # print(N:, N) omega = omega_plus[0:N, :] u_hat = np.zeros([T, K]) u_hat = u_hat + 0.j u_hat[T//2: T: 1, :] = np.squeeze(u_hat_plus[N-1, T//2: T: 1, :]) for i in range(0, T//2): u_hat[T // 2 - i, :] = np.squeeze(np.conj(u_hat_plus[N-1, T // 2 + i, :])) u_hat[0, :] = np.conj(u_hat[-1, :]) # print(u_hat[259, 2]) u = np.zeros([K, len(t)]) u = u + 0.j for k in range(0, K): u[k, :] = np.real(np.fft.ifft(np.fft.ifftshift(u_hat[:, k]))) # print(u[1, 99]) # 移除镜像部分 u = u[:, T//4:3*T//4] # print(u[1, 500:505]) # 重新计算频谱 u_hat = np.zeros([T//2, K]) u_hat = u_hat + 0.j for k in range(0, K): u_hat[:, k] = np.fft.fftshift(np.fft.fft(u[k, :])) # print(u_hat[201:207, 1],u_hat[199:206, 1]) return u, omega, u_hat
(系统自动生成,下载前可以参看下载内容)
下载文件列表
相关说明
- 本站资源为会员上传分享交流与学习,如有侵犯您的权益,请联系我们删除.
- 本站是交换下载平台，提供交流渠道，下载内容来自于网络，除下载问题外，其它问题请自行百度。
- 本站已设置防盗链，请勿用迅雷、QQ旋风等多线程下载软件下载资源，下载后用WinRAR最新版进行解压.
- 如果您发现内容无法下载，请稍后再次尝试；或者到消费记录里找到下载记录反馈给我们.
- 下载后发现下载的内容跟说明不相乎，请到消费记录里找到下载记录反馈给我们，经确认后退回积分.
- 如下载前有疑问，可以通过点击"提供者"的名字，查看对方的联系方式，联系对方咨询.
相关搜索:
