基于小波分析的GNSS信号周跳检测（Python）

摘要：import scipy.ioimport sysimport osimport matplotlib.pyplot as pltsys.path.append('../')from wavelet_1D import Wavelet_1Dplt.rc('te

import scipy.ioimport sysimport osimport matplotlib.pyplot as pltsys.path.append('../')from wavelet_1D import Wavelet_1Dplt.rc('text', usetex=True)plt.rc('font', family='serif')SAVE_OUNTPUTS = Falseif SAVE_OUNTPUTS:OUTPUT_FORMAT = '.pdf'OUTPUT_DIR = './Outputs'if not os.path.isdir(OUTPUT_DIR):os.mkdir(OUTPUT_DIR)mat_file = scipy.io.loadmat('CycleSlip2_L1_PRN03_GPS_512.mat')Signal = mat_file['Signal'] Time = mat_file['Time'] plt.figure(figsize=(6,4))# plt.plot(Time, Signal, color='k', linewidth=2.4)plt.plot(Time, Signal, color='k', marker='.', linestyle='', markersize=4)plt.gridplt.title('GPS PRN03')plt.xlabel('Time [hour of day]')plt.ylabel('L1 [meters]')plt.tight_layoutif SAVE_OUNTPUTS:plt.savefig(OUTPUT_DIR+'/CycleSlip_MainSignal'+OUTPUT_FORMAT, bbox_inches='tight')plt.show

S_haar, D_haar = Wavelet_1D(Signal, 'Haar', 3)S_db4, D_db4 = Wavelet_1D(Signal, 'Db4', 3)S_db6, D_db6 = Wavelet_1D(Signal, 'Db6', 3)S_mh, D_mh = Wavelet_1D(Signal, 'MexicanHat', 3)S_sym2, D_sym2 = Wavelet_1D(Signal, 'sym2', 3)import matplotlib.pyplot as pltimport numpy as npfrom scipy.signal import find_peaksn = len(Signal)F = np.zeros((n, n-5))g = np.array([[1], [-1], [-2], [2], [1], [-1]])for i in range(n-5):F[:, i] = np.hstack((np.zeros((1, i)) , g.T, np.zeros((1, n-i-6))))Y = F.T@Signalpeaks, _ = find_peaks(Y.flatten, height=500)plt.figure(figsize=(6,4))plt.plot(np.linspace(0, len(Y), len(Y)), Y, color='b')plt.plot(np.linspace(0, len(Y), len(Y))[peaks], Y[peaks], 'x', color='r', markersize=7) plt.gridplt.ylabel('Wavelet Transform')plt.xlabel('Time Index')plt.tight_layoutplt.showcs_idx_values = ', '.join(str(idx) for idx in peaks+2)time_values = ', '.join(f'{time[0]:.2f}' for time in Time[peaks+2]) print(f'Using Emperical Wavlet Transform - Cycle slip was detected in the time indices of {cs_idx_values} ({time_values} hour of day)')

import matplotlib.pyplot as pltplt.rc('text', usetex=True)plt.rc('font', family='serif')# %matplotlib tkn = len(S_haar)plt.figure(figsize=(12,9))for i in range(n):plt.subplot(n+1, 2, i*2+1)plt.plot(Time, Signal,linestyle='-', color='k', label = 'Original Signal')plt.plot(Time, S_haar[i],linestyle='-', color='r', label = f'Smoothed Signal ($S^{i+1}$)')plt.ylabel('GPS L1 [meters]')plt.gridplt.title(f'Smooth - Level {i+1} ($S^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]')else:plt.gca.set_xticklabels()plt.subplot(n+1, 2, i*2+2)plt.plot(Time, D_haar[i],linestyle='-', color='b', label = f'Details ($D^{i+1}$)')plt.gridplt.title(f'Detail - Level {i+1} ($D^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]') else:plt.gca.set_xticklabels()plt.suptitle('Haar Wavelet', fontsize=20)plt.tight_layoutif SAVE_OUNTPUTS:plt.savefig(OUTPUT_DIR+'/CycleSlip_Haar'+OUTPUT_FORMAT, bbox_inches='tight')plt.show

import matplotlib.pyplot as pltplt.rc('text', usetex=True)plt.rc('font', family='serif')# %matplotlib tkn = len(S_db4)plt.figure(figsize=(12,9))for i in range(n):plt.subplot(n+1, 2, i*2+1)plt.plot(Time, Signal,linestyle='-', color='k', label = 'Original Signal')plt.plot(Time, S_db4[i],linestyle='-', color='r', label = f'Smoothed Signal ($S^{i+1}$)')plt.ylabel('GPS L1 [meters]')plt.gridplt.title(f'Smooth - Level {i+1} ($S^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]')else:plt.gca.set_xticklabels()plt.subplot(n+1, 2, i*2+2)plt.plot(Time, D_db4[i],linestyle='-', color='b', label = f'Details ($D^{i+1}$)')plt.gridplt.title(f'Detail - Level {i+1} ($D^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]') else:plt.gca.set_xticklabels()plt.suptitle('Daubechies4 Wavelet', fontsize=20)plt.tight_layoutif SAVE_OUNTPUTS:plt.savefig(OUTPUT_DIR+'/CycleSlip_Daubechies4'+OUTPUT_FORMAT, bbox_inches='tight')plt.show

import matplotlib.pyplot as pltplt.rc('text', usetex=True)plt.rc('font', family='serif')# %matplotlib tkn = len(S_db6)plt.figure(figsize=(12,9))for i in range(n):plt.subplot(n+1, 2, i*2+1)plt.plot(Time, Signal,linestyle='-', color='k', label = 'Original Signal')plt.plot(Time, S_db6[i],linestyle='-', color='r', label = f'Smoothed Signal ($S^{i+1}$)')plt.ylabel('GPS L1 [meters]')plt.gridplt.title(f'Smooth - Level {i+1} ($S^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]')else:plt.gca.set_xticklabels()plt.subplot(n+1, 2, i*2+2)plt.plot(Time, D_db6[i],linestyle='-', color='b', label = f'Details ($D^{i+1}$)')plt.gridplt.title(f'Detail - Level {i+1} ($D^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]') else:plt.gca.set_xticklabels()plt.suptitle('Daubechies6 Wavelet', fontsize=20)plt.tight_layoutif SAVE_OUNTPUTS:plt.savefig(OUTPUT_DIR+'/CycleSlip_Daubechies6'+OUTPUT_FORMAT, bbox_inches='tight')plt.show

import matplotlib.pyplot as pltplt.rc('text', usetex=True)plt.rc('font', family='serif')# %matplotlib tkn = len(S_mh)plt.figure(figsize=(12,9))for i in range(n):plt.subplot(n+1, 2, i*2+1)plt.plot(Time, Signal,linestyle='-', color='k', label = 'Original Signal')plt.plot(Time, S_mh[i],linestyle='-', color='r', label = f'Smoothed Signal ($S^{i+1}$)')plt.ylabel('GPS L1 [meters]')plt.gridplt.title(f'Smooth - Level {i+1} ($S^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]')else:plt.gca.set_xticklabels()plt.subplot(n+1, 2, i*2+2)plt.plot(Time, D_mh[i],linestyle='-', color='b', label = f'Details ($D^{i+1}$)')plt.gridplt.title(f'Detail - Level {i+1} ($D^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]') else:plt.gca.set_xticklabels()plt.suptitle('Mexican Hat Wavelet', fontsize=20)plt.tight_layoutif SAVE_OUNTPUTS:plt.savefig(OUTPUT_DIR+'/CycleSlip_MexicanHat'+OUTPUT_FORMAT, bbox_inches='tight')plt.show

import matplotlib.pyplot as pltplt.rc('text', usetex=True)plt.rc('font', family='serif')# %matplotlib tkn = len(S_sym2)plt.figure(figsize=(12,9))for i in range(n):plt.subplot(n+1, 2, i*2+1)plt.plot(Time, Signal,linestyle='-', color='k', label = 'Original Signal')plt.plot(Time, S_sym2[i],linestyle='-', color='r', label = f'Smoothed Signal ($S^{i+1}$)')plt.ylabel('GPS L1 [meters]')plt.gridplt.title(f'Smooth - Level {i+1} ($S^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]')else:plt.gca.set_xticklabels()plt.subplot(n+1, 2, i*2+2)plt.plot(Time, D_sym2[i],linestyle='-', color='b', label = f'Details ($D^{i+1}$)')plt.gridplt.title(f'Detail - Level {i+1} ($D^{i+1}$)')if i == n-1:plt.xlabel('Time [hour of day]') else:plt.gca.set_xticklabels()plt.suptitle('Symlet2 Wavelet', fontsize=20)plt.tight_layoutif SAVE_OUNTPUTS:plt.savefig(OUTPUT_DIR+'/CycleSlip_Symlet2'+OUTPUT_FORMAT, bbox_inches='tight')plt.show

知乎学术咨询：

https://www.zhihu.com/consult/people/792359672131756032?isMe=1

担任《Mechanical System and Signal Processing》《中国电机工程学报》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测。

分割线分割线

集合高阶统计量和小波块阈值的非平稳信号降噪方法-以地震信号为例（MATLAB）

完整数据可通过知乎学术咨询获得

一种新的类谱峭度算法的旋转机械故障诊断模型（Python）

完整数据可通过知乎学术咨询获得

采用8种方法对一维信号进行降噪（Python）

NS: noisy signalS: original siganlmean filter: ws = window sizemedian filter:average filter: ns = number of noisy signal(different)bandpass filter: l = low cut-off frequency, h = high ...threshold filter: r = ratio(max abs(fft) / min ...)wavelet filter: a = thresholdstd filter:NN: neural network