iT邦幫忙

0

[Python] vibration analysis

分析如震動、音頻等資訊時,時間序列的資料往往分析不出太多資訊
以下以快速傅立葉轉換(Fast Fourier Transform, FFT)
可以繪製出
* 原始資訊波形
* RMS波形
* 時頻圖
* 時頻圖(log Data)

vibra_ana(10240,df_vib.DateTime,[df_vib.X,df_vib.Y,df_vib.Z])
def vibra_ana(Fs,t,x):
    
    print('from {} to {}'.format(t[0],t[len(t)-1]))
    
    l=len(x)
    figsize = (10 * l, 20 )
    row = 5
    col = l

    plt.figure(figsize=figsize)
    
 #Determine variables
    N = np.int(np.prod(t.shape))#length of the array
    T = 1/Fs
    T = (t[1]-t[0]).microseconds*0.000001#timedelta(seconds=0.2).total_seconds()
    Fs = 1/T 	#sample rate (Hz)
    print("Fs :{} , T :{}".format(Fs,T))
    print("# Samples:",N)
 
    for index_ele , ele in enumerate(x) :
        index = 0 * l + index_ele + 1 
        
        #Plot Data
        tic = time.clock()
        plt.subplot(row,col,index)   
        plt.plot(t, ele)
        plt.xlabel('Time (seconds)')
        plt.ylabel('Accel (g)')
        plt.title('RAW Data - ' )
        plt.grid()
        toc = time.clock()
        print("Plot Time:",toc-tic)
        
        index = 1 * l + index_ele + 1 

        #Compute RMS and Plot
        tic = time.clock()
        w = np.int(np.floor(Fs/10)); #width of the window for computing RMS
        print('w :{}'.format(w))
        steps = np.int_(np.floor(N/w)) #Number of steps for RMS
        t_RMS = np.zeros((steps,1)); #Create array for RMS time values
        x_RMS = np.zeros((steps,1)); #Create array for RMS values
        for i in range (0, steps):
#             t_RMS[i] = np.mean([(t[(i*w)]-t[0]).microseconds,(t[((i+1)*w-1)]-t[0]).microseconds])
            t_RMS[i] = np.mean([(t[(i*w)]-t[0]).seconds,(t[((i+1)*w-1)]-t[0]).seconds])
    #         t_RMS[i] = np.mean(t[(i*w):((i+1)*w)]);
            x_RMS[i] = np.sqrt(np.mean(ele[(i*w):((i+1)*w)]**2));  
        plt.subplot(row,col,index)   
        plt.plot(t_RMS, x_RMS)
        plt.xlabel('Time (seconds)')
        plt.ylabel('RMS Accel (g)')
        plt.title('RMS - ' )
        plt.grid()
        toc = time.clock()
        print("RMS Time:",toc-tic)
        index = 2 * l + index_ele + 1 


        #Compute and Plot FFT
        tic = time.clock()
        plt.subplot(row,col,index)   
        xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
        yf = fft(ele)
        plt.plot(xf, 2.0/N * np.abs(yf[0:np.int(N/2)]))
        plt.grid()
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Accel (g)')
        plt.title('FFT - ' )
        toc = time.clock()
        print("FFT Time:",toc-tic)
        index = 3 * l + index_ele + 1 


        
        # Compute and Plot Spectrogram
        tic = time.clock()
        plt.subplot(row,col,index)   
        f, t2, Sxx = signal.spectrogram(ele, Fs )#, nperseg = Fs/4)
        ax=plt.pcolormesh(t2, f, np.log(Sxx))
#         plt.colorbar(ax)
        plt.ylabel('Frequency [Hz]')
        plt.xlabel('Time [sec]')
        plt.title('Spectrogram(log) - ' )
        toc = time.clock()
        print("Spectrogram Time:",toc-tic)
        index = 4 * l + index_ele + 1 


    # Compute and Plot Spectrogram
        tic = time.clock()
        plt.subplot(row,col,index)   
        plt.pcolormesh(t2, f, (Sxx))
        plt.ylabel('Frequency [Hz]')
        plt.xlabel('Time [sec]')
        plt.title('Spectrogram - ' )
        toc = time.clock()
        print("Spectrogram Time:",toc-tic)

    
    plt.tight_layout()
    plt.show()

尚未有邦友留言

立即登入留言