一维信号的时频分析（Python）

代码较为简单，很容易读懂。

Importing the required libraries

import os
import numpy as np
import pywt
import pandas as pd
import pickle as pkl
from matplotlib import pyplot as plt

Parameters or Required Variables

DATA_POINTS_PER_FILE = 2560
TIME_PER_REC = 0.1
SAMPLING_FREQ = 25600 # 25.6 KHz
SAMPLING_PERIOD = 1.0/SAMPLING_FREQ


WIN_SIZE = 20
WAVELET_TYPE = 'morl'

Helper Functions

def load_df(pkz_file):
    with open(pkz_file, 'rb') as f:
        df=pkl.load(f)
    return df

# perform CWT on 1d signals and return 2d feature image
def extract_feature_image(ind, feature_name='horiz accel'):
    data_range = df_row_ind_to_data_range(ind)
    data = df[feature_name].values[data_range[0]:data_range[1]]
    # use window to process(= prepare, develop) 1d signal
    data = np.array([np.mean(data[i:i+WIN_SIZE]) for i in range(0, DATA_POINTS_PER_FILE, WIN_SIZE)])  
    # perform cwt on 1d data
    coef, _ = pywt.cwt(data, np.linspace(1,128,128), WAVELET_TYPE)  
    # transform to power and apply logarithm ?!
    coef = np.log2(coef**2+0.001) 
    # normalize coef
    coef = (coef - coef.min())/(coef.max() - coef.min()) 
    return coef

main_dir = ''

pkz_file=main_dir+'bearing1_3.pkz'
df=load_df(pkz_file)
df.head()

Printing total data points and total no. of data files present in a pickle(.pkz) file

no_of_rows = df.shape[0]
no_of_files = int(no_of_rows / DATA_POINTS_PER_FILE)
print(no_of_rows, no_of_files)

Plotting 1D vibration signals(both horiz accel and vert accel)

plt.plot(range(no_of_rows), df['horiz accel'])
plt.show()
plt.plot(range(no_of_rows), df['vert accel'], 'r')
plt.show()

signal processing = Extracting Time-Frequency Domain feature images

no_of_samples=5
fig, ax = plt.subplots(2, no_of_samples, figsize=[20,8])
ax[0,0].set_ylabel('horiz accel features image')
ax[1,0].set_ylabel('vert accel features image')
'''
dividing the feature images into 5 samples
'''
for i, p in enumerate(np.linspace(0,1,no_of_samples)):
    ind = int((no_of_files-1)*p)




  #extracting and plotting horizontal acceleration feature images (horiz accel feature images) for 5 samples (0.00, 0.25, 0.50, 0.75, 1.00).
  #horiz accel => horizontal acceleration vibration signal


    coef = extract_feature_image(ind, feature_name='horiz accel')
    ax[0,i].set_title('{0:.2f}'.format(p))
    im = ax[0,i].imshow(coef, cmap='coolwarm')
    fig.colorbar(im, ax=ax[0,i], fraction=0.046, pad=0.04)


  
  #extracting and plotting vertical acceleration feature images (vert accel feature images) for 5 samples (0.00, 0.25, 0.50, 0.75, 1.00). 
  #vert accel => vertical acceleration vibration signal
  
    coef = extract_feature_image(ind, feature_name='vert accel')
    ax[1,i].set_title('{0:.2f}'.format(p))
    im = ax[1,i].imshow(coef, cmap='coolwarm')
    fig.colorbar(im, ax=ax[1,i], fraction=0.046, pad=0.04)
  
'''
The tight_layout() function in pyplot module of matplotlib library is used to automatically adjust subplot parameters to give specified padding.
padding = adding space (adding required white space)
'''
plt.tight_layout()
'''
The show() function in pyplot module of matplotlib library is used to display all figures.
'''
plt.show()

工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。