Python环境下一种改进的基于梯度下降的自适应短时傅里叶变换

在数字信号处理技术中，傅里叶变换及其逆变换是一种信号时频分析方法。该方法将信号的时域描述及频域描述联系在一起，时域信号可通过正变换转变为频域信号，频域信号可通过逆变换转变为时域信号进行分析。但傅里叶变换及其逆变换是一种信号的整体变换，要么完全在时域内进行分析处理，要么完全在频域内进行分析处理，无法得到信号频谱含量随时间的变化规律，即传统傅里叶变换及其逆变换存在不适用于非平稳信号、没有局限性及时域与频域分割的缺陷。要克服上述缺陷，就需要寻找联合时频方法。

1946年，GaborD.提出了短时傅里叶变换，该变换可将信号的时域分析和频域分析结合起来，其变换结果既反映了信号的频率内容，又反映了频率内容随时间变化的规律。短时傅里叶变换的过程为：用一个时宽足够窄的窗函数与时间信号相乘，窗内的信号近似视为平稳信号，然后在窗内进行傅里叶变换，得到信号的瞬时频谱。随着截取窗在时间轴上移动，得到整个时间域上的频谱。

短时傅里叶变换是一种线性的联合时频分析方法，其采用固定的窗函数对信号进行加窗处理，该变换克服了传统傅里叶变换的缺陷，同时易于实现。与傅里叶分析不同，短时傅里叶变换中的窗函数种类繁多，为将短时傅里叶变换的效果最优化，窗函数的选取显得尤为重要。实际应用的基本窗函数，主要分为以下类型：应用幂次时间函数的幂窗、应用三角时间函数的三角函数窗及应用指数时间函数的指数窗。

短时傅里叶变换中，采用窗函数对信号进行截取会造成频谱泄露，即某一频率的信号能量扩散到其相邻频率点。频谱泄漏与窗函数频谱的两侧旁瓣有关，从保持信号最大信息和消除旁瓣的综合效果出发，窗函数选取的基本原则为：（1）尽可能使窗函数频谱中的主瓣宽度足够窄，以获得较陡的过渡带，能量尽可能集中在主瓣内，提高频谱分析时的频率分辨率；（2）旁瓣高度尽量小且随频率尽快衰减，以提高阻带的衰减，减小频谱分析时的泄露失真。但主瓣既窄，旁瓣又小衰减又快的窗函数很难找到，所以需要根据信号的性质和研究目的来选取合适的窗函数。

为了自适应的对窗长进行调整，提出一种基于梯度下降的自适应短时傅里叶变换，主要贡献是窗口长度优化的新范例，包括修改STFT算子的定义，使窗长成为一个连续的参数，进而使得谱图可以被微分，运行环境为Python，采用Jupyter Notebook文件格式。部分模块如下：

from math import pi
import matplotlib.pyplot as plt
import numpy as np
import torch
import torchaudio

部分代码如下：

sr = 1e4
t = torch.arange(0, 2.5, 1/sr)
f = torch.sin(2*pi*t) * 1e2 + 1e2 * torch.ones_like(t) + 5e1 * t
x = (torch.sin(torch.cumsum(f, dim=0) / 2e2) + 0.1 *torch.randn(t.shape))[None, :]
x += torch.sin(torch.cumsum(1e2*5 * torch.ones_like(t), dim=0) / 2e2)
x = x.to(device)
print(x.shape)

plt.plot(f)

spec = torchaudio.transforms.Spectrogram(n_fft=1_000, win_length=1_000, hop_length=100)(x.cpu())
plt.imshow(spec[0].log().cpu(), aspect='auto', origin='lower', cmap='jet', extent=[0,spec.shape[-1], 0, spec.shape[-2]])
plt.ylabel('frequencies (Hz)', fontsize=18)
plt.xlabel('frames', fontsize=18)

出图如下：