第二章：整車發動機激勵--快速傅里葉變換

摘要：第一章介紹了缸壓曲線和曲軸轉速轉化成曲軸中心點的時間-集中力；但是這個集中力並不是我們可以直接使用的頻域-集中力；需要經過快速傅里葉變換，將時域力轉化成頻域力。

• 快速傅里葉變換的理解

　　　　還是老規矩，我們需要什么？我們需要從時域到頻域的轉換。怎么理解呢？當然是先找一個已知的頻率來看看是什么樣的。從網上盜了一張動態圖，從圖中我們已知各個三角函數的周期（頻率），可以做出復雜的曲線。【先看成果，我們的目的是從復雜的曲線中分離出各個周期（頻率）】

"采樣點選取,0-2間取2^n個采樣點,時域橫坐標，單位s"
n=10
x=np.linspace(0,2,2**n)
# print(x)
"構造信號曲線，其中包含頻率1Hz，2Hz，3Hz，4Hz"
y=2*np.sin(2*np.pi*1*x)/np.pi-2*np.sin(2*np.pi*2*x)/(2*np.pi)+2*np.sin(2*np.pi*3*x)/(3*np.pi)\
  -2*np.sin(2*np.pi*4*x)/(4*np.pi)

plt.plot(x, y)
plt.show()

 

• 快速傅里葉變換

 

import numpy as np
from scipy.fftpack import  fft,ifft
import matplotlib.pyplot as plt
import seaborn

"采樣點選取,0-2間取2^n個采樣點,時域橫坐標，單位s"
n=10
x=np.linspace(0,2,2**n)
# print(x)
"構造信號曲線，其中包含頻率1Hz，2Hz，3Hz，4Hz"
y=2*np.sin(2*np.pi*1*x)/np.pi-2*np.sin(2*np.pi*2*x)/(2*np.pi)+2*np.sin(2*np.pi*3*x)/(3*np.pi)\
  -2*np.sin(2*np.pi*4*x)/(4*np.pi)

plt.plot(x, y)
plt.show()
yy=fft(y) #快速傅里葉變換
yreal = yy.real  #獲取實部
ymiag = yy.imag  #獲取虛部
yf=abs(fft(y)) #取絕對值
yf1 = abs(fft(y))/len(x) #歸一化處理
yf2 = yf1[range(int(len(x)/2))] ##對稱性只取一半區間

"采樣頻率"
dt=x[1]-x[0]
fs = 1/dt
freqs = fs/2*np.linspace(0,1,int(2**n/2))
print(len(yf2))
print(len(freqs))
plt.subplot(211)
plt.plot(x, y)
plt.title('Original wave')
plt.subplot(212)
plt.plot(freqs[0:10], yf2[0:10])
plt.title('FFT wave')
plt.show()

從及結果可以看出，通過快速傅里葉變換，找到了四個階次的頻率，分別是1Hz，2Hz，3Hz，4Hz，與構造的信號函數頻率一一對應。

 

• 傅里葉變換處理發動機力（時域-頻域）

　　　　上述闡述已知信號源的傅里葉變換，其實信號源是否已知並不影響傅里葉變換的結果，構造已知信號的函數，僅僅是為了直觀的理解這個變換輸入輸出。【這里並未對傅里葉變換的底層實現進行闡述，僅僅說明如何在發動機載荷分解過程中需要用到的功能】。

　　　　回到上一章接種獲取到的發動機曲軸中心位置的時間-力的數據。每一組時間-力（力矩）都如同上訴的信號源（只是未知信號）。

　　　　

 

 　　　　

 

 

 

import numpy as np
from scipy.fftpack import  fft,ifft
import matplotlib.pyplot as plt
import seaborn
from scipy.interpolate import interp1d
"讀取時間載荷"
data = np.loadtxt(r'TimeLoad\A750 _transpose.csv',dtype=np.float, delimiter=",")
"獲取時間采樣點"
t=data[0]
"獲取z向時間載荷"
z=data[3]
print(t)
print(z)
"原始數據采樣，采樣點數量按照2**n個進行，運用插值法，獲取格則采樣數據,"
n=13
NFFT = 2**n
x, delta_step = np.linspace(t[0], t[len(t) - 1], NFFT, endpoint=True, retstep=True)
y_curve = interp1d(t,z,kind='linear')
y = y_curve(x)

plt.plot(x, y)
plt.show()
yy=fft(y) #快速傅里葉變換
yreal = yy.real  #獲取實部
ymiag = yy.imag  #獲取虛部
yf=abs(fft(y)) #取絕對值
yf1 = abs(fft(y))/len(x) #歸一化處理
yf2 = yf1[range(int(len(x)/2))] ##對稱性只取一半區間

"采樣頻率"
dt=x[1]-x[0]
fs = 1/dt
freqs = fs/2*np.linspace(0,1,int(2**n/2))
print(len(yf2))
print(len(freqs))
plt.subplot(211)
plt.plot(x, y)
plt.title('Original wave')
plt.subplot(212)
plt.plot(freqs[0:100], yf2[0:100])
plt.title('FFT wave')
plt.show()

 

在怠速750r/min時，在25Hz的激勵是62.157N。