四虎免费看片,一级理伦,亚洲综合涩

1 短時傅里葉變換STFT原理介紹

1.1 傅里葉變換的本質

傅里葉變換的基本思想：

將信號分解成一系列不同頻率的連續正弦波的疊加；
或者說，將信號從時間域轉換到頻率域。

由于傅里葉變換是對整個信號進行變換，將整個信號從時域轉換到頻域，得到一個整體的頻譜；丟掉了時間信息，無法根據傅立葉變換的結果判斷一個特定信號在什么時候發生；所以傅里葉變換缺乏時頻分析能力、多分辨率分析能力，難以分析非平穩信號。

1.2 STFT概述

短時傅里葉變換（Short-Time Fourier Transform，STFT）是一種將信號分解為時域和頻域信息的時頻分析方法。它通過將信號分成短時段，并在每個短時段上應用傅里葉變換來捕捉信號的瞬時頻率。即采用中心位于時間α的時間窗g(t-α)在時域信號上滑動，在時間窗g(t-α)限定的范圍內進行傅里葉變換，這樣就使短時傅里葉變換具有了時間和頻率的局部化能力，兼顧了時間和頻率的分析[1]。

使用窄窗，時間分辨率提高而頻率分辨率降低；
使用寬窗，頻率分辨率提高而時間分辨率降低。

比如用利用高斯窗STFT對非平穩信號進行分析：

1.3 STFT的原理和過程

1.3.1 時間分割

在短時傅里葉變換（STFT）中，首先將信號分割成短時段。這個過程通常通過使用窗口函數來實現，窗口函數是一個在有限時間內非零，而在其他時間上為零的函數。常見的窗口函數有矩形窗、漢明窗、漢寧窗等。通過將窗口函數應用于信號，可以將信號分成許多短時段。

1.3.2 傅里葉變換

對于每個短時段，都會進行傅里葉變換。傅里葉變換是一種將信號從時域（時間域）轉換為頻域（頻率域）的方法。在這個上下文中，它用于分析每個短時段內信號的頻率成分。傅里葉變換將信號表示為不同頻率的正弦和余弦函數的組合。

1.3.3 時頻圖：

將每個短時段的傅里葉變換結果排列成一個矩陣，構成了時頻圖。時頻圖的橫軸表示時間，縱軸表示頻率，而每個點的強度表示對應頻率在對應時刻的幅度。時頻圖提供了一種直觀的方式來觀察信號在時間和頻率上的變化。

1.4 公式表示

STFT的數學表達式如下：

其中，

2 基于Python的STFT實現與參數對比

在 Python 中，使用 scipy 庫來實現短時傅里葉變換（STFT）

2.1 代碼示例

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import stft
# 生成三個不同頻率成分的信號
fs = 1000  # 采樣率
t = np.linspace(0, 1, fs, endpoint=False)  # 時間


# 第一個頻率成分
signal1 = np.sin(2 * np.pi * 50 * t)
# 第二個頻率成分
signal2 = np.sin(2 * np.pi * 120 * t)
# 第三個頻率成分
signal3 = np.sin(2 * np.pi * 300 * t)


# 合并三個信號
signal = np.concatenate((signal1, signal2, signal3))


# 進行短時傅里葉變換  
f, t, spectrum = stft(signal, fs, nperseg=100, noverlap=50)
# 繪制時頻圖
plt.pcolormesh(t, f, np.abs(spectrum), shading='gouraud')
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

參數解釋

nperseg：表示每個段的長度，即窗口大小
noverlap：表示相鄰段的重疊部分
f：是頻率數組，表示頻譜的頻率分量
t：是時間數組，表示每個時間段的起始時間
spectrum：是頻譜矩陣，spectrum[f, t] 表示在頻率 f 處于時間 t 時的頻譜強度

可以使用 np.abs(spectrum) 獲取頻譜的幅度，表示在不同頻率和時間上的信號強度

2.2 參數選擇和對比

2.2.1 nperseg（窗口長度）：

nperseg 定義了每個時間段內的信號長度，也就是說，它規定了窗口的大小。
**較小的 **nperseg 可以提供更高的時間分辨率，因為窗口更短，可以捕捉到信號的更快變化。但頻率分辨率較差。
較大的 nperseg 提供更高的頻率分辨率，但可能失去對信號快速變化的捕捉。
選擇適當的窗口長度要根據信號特性，如頻率變化的速率和信號長度等。

2.2.2 noverlap（重疊長度）：

noverlap 定義了相鄰時間段之間的重疊部分。
**較大的 **noverlap 可以提高時間分辨率，因為相鄰時間段之間有更多的共享信息。但可能導致頻譜圖中的泄漏（leakage）。
**較小的 **noverlap 可以提高頻率分辨率，減少泄漏，但時間分辨率可能下降。

2.2.3 選擇策略：

參數的選擇需要考慮到對信號分析的具體需求，平衡時間和頻率分辨率，嘗試不同的 nperseg 和 noverlap 值，觀察結果的變化。

對于 nperseg，選擇較小的值可能需要根據信號的特定頻率內容進行調整。確保窗口足夠短，以捕捉到頻率變化。
對于 noverlap，通常選擇 50% 的重疊是一個常見的起點

2.3 凱斯西儲大學軸承數據的加載

選擇正常信號和 0.021英寸內圈、滾珠、外圈故障信號數據來做對比

第一步，導入包，讀取數據

import numpy as np
from scipy.io import loadmat
import numpy as np
from scipy.signal import stft
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rc("font", family='Microsoft YaHei')


# 讀取MAT文件
data1 = loadmat('0_0.mat')  # 正常信號
data2 = loadmat('21_1.mat') # 0.021英寸 內圈
data3 = loadmat('21_2.mat') # 0.021英寸 滾珠
data4 = loadmat('21_3.mat') # 0.021英寸 外圈
# 注意，讀取出來的data是字典格式，可以通過函數type(data)查看。

第二步，數據集中統一讀取驅動端加速度數據，取一個長度為1024的信號進行后續觀察和實驗

# DE - drive end accelerometer data 驅動端加速度數據
data_list1 = data1['X097_DE_time'].reshape(-1)
data_list2 = data2['X209_DE_time'].reshape(-1)  
data_list3 = data3['X222_DE_time'].reshape(-1)
data_list4 = data4['X234_DE_time'].reshape(-1)
# 劃窗取值（大多數窗口大小為1024）
data_list1 = data_list1[0:1024]
data_list2 = data_list2[0:1024]
data_list3 = data_list3[0:1024]
data_list4 = data_list4[0:1024]

第三步，進行數據可視化

plt.figure(figsize=(20,10))


plt.subplot(2,2,1)
plt.plot(data_list1)
plt.title('正常')
plt.subplot(2,2,2)
plt.plot(data_list2)
plt.title('內圈')
plt.subplot(2,2,3)
plt.plot(data_list3)
plt.title('滾珠')
plt.subplot(2,2,4)
plt.plot(data_list4)
plt.title('外圈')


plt.show()

2.4 STFT與參數選擇

2.4.1 基于重疊比例為0.5，選擇內圈數據比較 STFT 的不同尺度：16、32 、64、128

from scipy.signal import stft


# 設置STFT參數
window_size = 32  # 窗口大小
overlap = 0.5  # 重疊比例
# 計算重疊的樣本數
overlap_samples = int(window_size * overlap)
frequencies1, times1, magnitude1 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples)


# 設置STFT參數
window_size = 64  # 窗口大小
overlap = 0.5  # 重疊比例
# 計算重疊的樣本數
overlap_samples = int(window_size * overlap)
frequencies2, times2, magnitude2 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples)


# 設置STFT參數
window_size = 128  # 窗口大小
overlap = 0.5  # 重疊比例
# 計算重疊的樣本數
overlap_samples = int(window_size * overlap)
frequencies3, times3, magnitude3 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples)


# 設置STFT參數
window_size = 256  # 窗口大小
overlap = 0.5  # 重疊比例
# 計算重疊的樣本數
overlap_samples = int(window_size * overlap)
frequencies4, times4, magnitude4 = stft(data_list2, nperseg=window_size, noverlap=overlap_samples

數據可視化

plt.figure(figsize=(20,10), dpi=100)
plt.subplot(2,2,1)
plt.pcolormesh(times1, frequencies1, np.abs(magnitude1), shading='gouraud')
plt.title('尺度 16-內圈')
plt.subplot(2,2,2)
plt.pcolormesh(times2, frequencies2, np.abs(magnitude2), shading='gouraud')
plt.title('尺度 32-內圈')
plt.subplot(2,2,3)
plt.pcolormesh(times3, frequencies3, np.abs(magnitude3), shading='gouraud')
plt.title('尺度 64-內圈')
plt.subplot(2,2,4)
plt.pcolormesh(times4, frequencies4, np.abs(magnitude4), shading='gouraud')
plt.title('尺度 128-內圈')
plt.show()