import cmathimport matplotlib.pyplot as pltimport numpy as npdef create_signal(frequency, time):sin = np.sin(2 * np.pi * (frequency * time))sin2 = np.sin(2 * np.pi * (2 * frequency * time))sin3 = np.sin(2 * np.pi * (3 * frequency * time))return sin + sin2 + sin3def calculate_centre_of_gravity(mult_signal):x_centre = np.mean([x.real for x in mult_signal])y_centre = np.mean([x.imag for x in mult_signal])return x_centre, y_centredef calculate_sum(mult_signal):x_sum = np.sum([x.real for x in mult_signal])y_sum = np.sum([x.imag for x in mult_signal])return x_sum, y_sumdef create_pure_tone(frequency, time):angle = -2 * np.pi * frequency * timereturn np.cos(angle) + 1j * np.sin(angle)def plot_fourier_transform(pure_tone_frequency, signal_frequency, time, plot_centre_of_gravity=False,plot_sum=False):# create sinusoid and signalpure_tone = create_pure_tone(pure_tone_frequency, time)signal = create_signal(signal_frequency, time)# multiply pure tone and signalmult_signal = pure_tone * signalX = [x.real for x in mult_signal]Y = [x.imag for x in mult_signal]plt.figure(figsize=(15, 10))plt.plot(X, Y, 'o')# calculate and plot centre of gravityif plot_centre_of_gravity:centre_of_gravity = calculate_centre_of_gravity(mult_signal)plt.plot([centre_of_gravity[0]], [centre_of_gravity[1]], marker='o', markersize=10, color="red")# calculate and plot sum if plot_sum:integral = calculate_sum(mult_signal)plt.plot([integral[0]], [integral[1]], marker='o', markersize=10, color="green")# set origin axesax = plt.gcaax.grid(True)ax.spines['left'].set_position('zero')ax.spines['right'].set_color('none')ax.spines['bottom'].set_position('zero')ax.spines['top'].set_color('none')if not plot_sum:plt.xlim(-3, 3)plt.ylim(-3, 3)plt.showdef plot_signal(signal, time):plt.figure(figsize=(15, 10))plt.plot(signal, time)plt.xlabel("Time")plt.ylabel("Intensity")plt.showtime = np.linspace(0, 10, 10000)signal = create_signal(frequency=1, time=time)plot_signal(time, signal)摘要：import cmathimport matplotlib.pyplot as pltimport numpy as npdef create_signal(frequency, time):sin = np.sin(2 * np.pi * (frequenc
知乎学术咨询：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1担任《Mechanical System and Signal Processing》《中国电机工程学报》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测。
分割线分割线
基于小波分析的时间序列降噪（Python,ipynb文件）
import osimport numpy as npfrom pywt import cwt, wavedec, swt, WaveletPacketfrom scipy.fft import rfft, irfft, rfftfreqimport matplotlib.pyplot as pltimport matplotlib.gridspec as gridspecplt.rcParams.update({"font.size": 12, "font.weight": 'bold', "lines.linewidth": 1.5})#Load the datapds = np.loadtxt(r"Data\pds_sim.txt")[:832, :]pds30 = np.loadtxt(r"Data\pds_sim30.txt")[:832]# Calculate number of sampling pint, sampling frequencyN, deltaT = len(pds[:, 0]), (pds[1, 0]-pds[0, 0])# Take Fourier transform of the datafft_pds, fft_pds30 = rfft(pds[:, 1]), rfft(pds30)fq = rfftfreq(N, deltaT)# Take Short Time Fourier Transform of the data# stft_fq, stft_t, stft_esr = stft(x=-1*data_esr[:, 1], fs=deltaT, window='hann', nperseg=256, noverlap=0)# _, _, stft_esr30 = stft(x=data_esr30, fs=deltaT, window='hann', nperseg=256, noverlap=0)seg = [0, 208, 416, 624]stft_pds, stft_pds30 = np.empty((len(seg), len(fft_pds)), dtype=np.complex64), np.empty((len(seg), len(fft_pds)), dtype=np.complex64)for i in range(len(seg)):tempd1, tempd2 = np.zeros((N)), np.zeros((N))tempd1[:208], tempd2[:208] = pds[seg[i]:seg[i]+208, 1], pds30[seg[i]:seg[i]+208]stft_pds[i, :], stft_pds30[i, :] = rfft(tempd1), rfft(tempd2)# Frequency filtering# copy the datafftbp_pds, fftbp_pds30 = np.copy(fft_pds), np.copy(fft_pds30)fftgp_pds, fftgp_pds30 = np.copy(fft_pds), np.copy(fft_pds30)# Perform low and high frequency filtering#fftbp_esr[(freq_esr (np.max(fq))/2)] = 0# generate a Gaussian bandpass functiongaussp = np.exp(-((fq - min(fq))/(2*16))**2) + np.exp(-((fq - max(fq))/(2*16))**2) # Possitive frequency#gaussn = np.exp(-((fq + min(fq))/(2*0.1))**2) + np.exp(-((fq + max(fq))/(2*0.1))**2) # Negative frequencygaussf = gaussp #+ gaussn # Only have possitive frequencies hence need to filter them onlyfftgp_pds30 = fftgp_pds30*gaussf# Reconstruct the data using inverse Fourier transformdatabp_pds, databp_pds30 = irfft(fftbp_pds), irfft(fftbp_pds30)datagp_pds, datagp_pds30 = irfft(fftgp_pds), irfft(fftgp_pds30)# Take FT of the denoised dataftbp_pds30, ftgp_pds30 = rfft(databp_pds30), rfft(datagp_pds30)# Tiime domain pds Plot and save the dataplt.figure(figsize=(4.5, 3), dpi=100, layout='constrained')# plt.plot(pds[:, 0], pds30[::-1], '-k', label='Noisy')plt.plot(pds[:, 0], pds30, '-r', label='Noisy')plt.plot(pds[:, 0], pds[:, 1], '-b', label='Noise-free')plt.legend(frameon=False, fontsize=10, loc='upper right')plt.xlim([(pds[0, 0]), (pds[831, 0])])plt.xlabel('Time ($\\mu s$)', fontsize=14, fontweight='bold'), plt.ylabel('Magnitude', fontsize=14, fontweight='bold')# plt.vlines(x=pds[207, 0], ymin=min(pds[:, 1])-0.01, ymax=max(pds[:, 1])+0.01, color='k', ls='--', lw=2)# plt.vlines(x=pds[415, 0], ymin=min(pds[:, 1])-0.01, ymax=max(pds[:, 1])+0.01, color='k', ls='--', lw=2)# plt.vlines(x=pds[623, 0], ymin=min(pds[:, 1])-0.01, ymax=max(pds[:, 1])+0.01, color='k', ls='--', lw=2)# plt.savefig('pds.png', dpi=400, bbox_inches='tight', pad_inches=0.02)plt.show完整代码：
信号降噪与时频分析（Python,ipynb文件）
完整代码：
海面水温序列的连续小波变换（Python，ipynb文件）
Plot the signal and its wavelet spectrum:
variance = np.std(s3, ddof=1) ** 2global_ws = variance * (np.sum(power, axis=1) / n) # time-averaged global spectrumfrom matplotlib import gridspecgs = gridspec.GridSpec(2, 2, width_ratios=[3,1], height_ratios=[1,3])#--- signalxlim = ([0, 130])plt.figure(figsize=(12, 10))plt.subplot(gs[0,0])plt.plot(t, s3)plt.xlim(xlim[:])plt.xlabel('time (year)')plt.ylabel('amplitudes')plt.title('3 sines, periods 1, 4, 16 years')#--- wavelet power spectrumplt3 = plt.subplot(gs[1,0])levels = np.array([2**i for i in range(-18,10)])CS = plt.contourf(t, period, np.log2(power), len(levels)) #*** or use 'contour'im = plt.contourf(CS, levels=np.log2(levels))plt.xlabel('time (year)')plt.ylabel('period (year)')plt.title('3 sines Morlet wavelet power spectrum')plt.xlim(xlim[:])# format y-scaleplt3.set_yscale('log', base=2, subs=None)plt.ylim([np.min(period), np.max(period)])ax = plt.gca.yaxisax.set_major_formatter(matplotlib.ticker.ScalarFormatter)plt3.ticklabel_format(axis='y', style='plain')plt3.invert_yaxis#--- global wavelet spectrumplt4 = plt.subplot(gs[1,1])plt.plot(global_ws, period)plt.xlabel('ampl. square')#plt.ylabel('period (year)')plt.title('global wavelet spectr.')plt.xlim([0, 1.25 * np.max(global_ws)])# format y-scaleplt4.set_yscale('log', base=2, subs=None)plt.ylim([np.min(period), np.max(period)])ax = plt.gca.yaxisax.set_major_formatter(matplotlib.ticker.ScalarFormatter)plt4.ticklabel_format(axis='y', style='plain')plt4.invert_yaxisplt.tight_layoutplt.show基于连续小波变换的信号滤波方法（Python，ipynb文件）
plt.figure(figsize=(12, 8))levels = np.array([2**i for i in range(-15,5)])with np.errstate(divide='ignore'):CS= plt.contourf(tint, period, np.log2(powerthr), len(levels)) im = plt.contourf(CS, levels=np.log2(levels))plt.xlabel('time (s)')plt.ylabel('period (s)')plt.title('tangential acceleration, filtered power spectrum')plt.xlim(xlim[:])plt.yscale('log', base=2, subs=None)plt.ylim([np.min(period), np.max(period)])ax = plt.gca.yaxisax.set_major_formatter(matplotlib.ticker.ScalarFormatter)plt.ticklabel_format(axis='y', style='plain')plt.gca.invert_yaxisplt.showplt.figure(figsize=(12, 8))plt.plot(tint, axr)plt.xlim(xlim[:])plt.xlabel('time (s)')plt.ylabel('a_x (m/s^2)')plt.title('tangential acceleration, inverse WT filtered')plt.showPeriodic acceleration transients at 33 s and 36 s are probably due to the slipping driven wheel. This is confirmed by zooming at these places:
plt.figure(figsize=(12, 6))plt.subplot(121)plt.plot(tint, axr)plt.xlim([33.2,33.7])plt.xlabel('time (s)')plt.ylabel('a_x (m/s^2)')plt.title('1st transient')plt.subplot(122)plt.plot(tint, axr)plt.xlim([35.8,36.3])plt.xlabel('time (s)')plt.ylabel('a_x (m/s^2)')plt.title('2nd transient')plt.show完整代码：
来源：完美教育
免责声明：本站系转载，并不代表本网赞同其观点和对其真实性负责。如涉及作品内容、版权和其它问题，请在30日内与本站联系，我们将在第一时间删除内容!
