音频处理5_STFT_窗口函数

短时傅里叶变换 (Short-Time Fourier Transform, STFT) 是一种时频谱转换算法，它通过在时间上移动窗口函数并计算窗口内信号的频谱来获得信号在时间和频率上的信息。填充信号可以确保每个窗口都有足够的数据进行频谱计算，特别是在窗口函数的边缘。

窗口函数主要用于信号处理中的短时傅里叶变换（STFT）、滤波器设计和其他需要对信号进行窗函数处理的场景， 可减少频谱泄漏，避免傅里叶变换中的频谱混叠

1.窗口函数

让窗口内的信号序列 乘以 一个系数系列，以调整信号的取值范围，
这个系数系列来自 窗口函数

2.频谱泄漏（Spectral Leakage）

频谱泄漏是信号的能量未集中在原始频率位置，而是扩散到了其他频率位置的现象

傅里叶变换是分段的，假设每段有一个主频率信号，且信号是周期的。但时间上截断信号会影响信号原有的周期性，这种非周期性会导致频域中主频能量扩散到其他频率分量。

频谱泄漏表现在频谱图中是 信号的主要频率成分（主瓣）周围会出现一些不希望出现的频率成分（旁瓣）。这些旁瓣是由于截断或不适当的窗口函数引起的，会影响频谱的清晰度和准确性。

窗口函数可以减少这种现象，通过对窗口内信号进行系数调整，窗口两端逐渐减小到零，从而平滑过渡，减少非周期性误差，即降低频率域中的副瓣（sidelobes）。

3.频谱特征

3.1 主瓣（mainlobe）

频谱中最大的峰值区域，它代表了信号的主要频率成分,

在应用窗口函数后的频谱图中，主瓣通常位于频谱的中心位置，并包含了信号的主要能量。

主瓣宽度决定了频谱的分辨率。主瓣越窄，频谱分辨率越高，可以更好地区分相近的频率成分。

3.2 旁瓣（Sidelobe）

旁瓣是频谱中位于主瓣两侧的较小峰值区域。旁瓣表示频谱泄漏，即信号的能量扩散到主频率之外的频率分量。

旁瓣越高，频谱泄漏越严重，这会影响频谱分析的精度。

旁瓣高度决定了频谱泄漏的程度。旁瓣越低，频谱泄漏越小，频谱分析的准确性越高。

4. 填充

截断前的操作，把周期函数的 “补全” or “延续”

填充是在窗口信号的两端添加一些额外的数据，以便在进行窗口函数处理时，避免边界处的窗口截断问题 ，截断会导致边界处的能量损失和频谱伪像。

下面是将[0.1, 0.2, 0.3, 0.4, 0.5] 填充为 [0.4, 0.3, 0.2, 0.1, 0.2, 0.3, 0.4, 0.5, 0.4, 0.3, 0.2]的代码：

signal = torch.tensor([0.1, 0.2, 0.3, 0.4, 0.5]) # 模拟信号

# STFT 参数
n_fft = 8
hop_size = 2

# 对信号进行填充
padded_signal = torch.nn.functional.pad(signal.unsqueeze(0), (int((n_fft-hop_size)/2), int((n_fft-hop_size)/2)), mode='reflect')
padded_signal = padded_signal.squeeze(0)

print("原始信号长度:", len(signal)) # 5
print("填充后的信号长度:", len(padded_signal)) # 11
print(padded_signal) 
# tensor([0.4000, 0.3000, 0.2000, 0.1000, 0.2000, 0.3000, 0.4000, 0.5000, 0.4000, 0.3000, 0.2000])

5. 汉宁函数 （Hann Window）

5.1 公式

$w(n)=0.5(1-cos\frac{N-1}{2\pi n})$

• n 是窗口内的样本索引, $0\le n\le N-1$

• N是总长度,

• w(n)是窗口函数的值, 它是一个序列，长度为N

5.2 函数特性

• 中心对称性: $w(n)\bar{w}(N-1-n)$

• 窗口的两端值为零，这有助于减小边界效应(平滑过渡，增加连续性)

• 窗口函数在中间达到最大值

6. 代码实例：

6.1 代码1

对原始5Hz正弦波乘以汉宁窗口系数，降低旁瓣（sidelobe），拉窄主瓣（mainlobe）：

结果如图：

旁瓣 1.0 => 0.4

主瓣 (50, 250) => (100, 200)

6.1 代码2

对一个440Hz信号转换为视频谱

• 没有padding和window的结果：
  

• 有padding，无window的结果：

效果好一些，没那么多隔断部分

• 有padding，有window的结果：

没有断层，特征反馈清晰

另外，傅里叶转换后结果是个复数，这个是复数实部的图:

这是复数虚部的：

7. 总结和代码

频谱特性：汉宁窗口的频谱具有较低的旁瓣（sidelobe），这意味着它能够有效地抑制频谱泄漏，主瓣（mainlobe）较宽，这会降低频率分辨率。

短时傅里叶变换：STFT通过对信号进行分段，对每个片段应用汉宁窗口，能减少频谱泄漏，获得更准确的频谱信息。

代码1：

import numpy as np
import matplotlib.pyplot as plt
import torch

# 生成一个示例信号
sampling_rate = 1000
t = np.linspace(0, 1, sampling_rate, endpoint=False)
freq = 5  # 频率为5Hz的正弦波
signal = np.sin(2 * np.pi * freq * t)

# 定义窗口函数
window_size = 256
hann_window = torch.hann_window(window_size) # $w(n)=0.5(1−cos\frac{N−1}{2πn}), n \in [0, N-1]$


# 应用窗口函数
windowed_signal = signal[:window_size] * hann_window.numpy()

# 绘制信号和窗口函数
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(signal[:window_size], label="Original Signal")
plt.plot(hann_window.numpy(), label="Hanning Window")
plt.legend()
plt.title("Original Signal and Hanning Window")

plt.subplot(2, 1, 2)
plt.plot(windowed_signal, label="Windowed Signal")
plt.legend()
plt.title("Windowed Signal")

plt.tight_layout()
plt.show()

代码2:

# STFT 参数
n_fft = 2048
hop_size = 512
win_size = 2048 # 窗口长度

# 模拟音频信号
sampling_rate = 16000
t = np.linspace(0, 1, sampling_rate)
freq = 440  # 频率为440Hz的正弦波
y = torch.tensor(np.sin(2 * np.pi * freq * t)).float()

# 生成汉宁窗口
hann_window = torch.hann_window(win_size)

# 对信号进行填充
y = torch.nn.functional.pad(y.unsqueeze(0), (int((n_fft-hop_size)/2), int((n_fft-hop_size)/2)), mode='reflect') # [17536]
y = y.squeeze(0)
print(y.shape)

# ------------------计算 STFT-1
spec = torch.stft(y, n_fft, hop_length=hop_size, win_length=win_size, 
    window=hann_window, center=True, pad_mode='reflect', normalized=False, 
    onesided=True, return_complex=True)
print(spec.shape)  # torch.Size([1025, 35]), 查看STFT结果的形状

# 转换为频谱幅度
#real = spec.real # 实部
#imag = spec.imag # 虚部
#spec_magnitude = torch.sqrt(real**2 + imag**2).numpy()
spec_magnitude = spec.abs()
print(spec)

#绘制频谱图
plt.figure(figsize=(10, 6))
plt.imshow(20 * np.log10(spec_magnitude + 1e-6), aspect='auto', origin='lower', extent=[0, t[-1], 0, sampling_rate / 2])
plt.colorbar(format='%+2.0f dB')
plt.title("STFT Magnitude Spectrogram")
plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.show()