傅里叶变换@（stft和istft）

Posted 2022-04-22 月疯

一、窗函数之短时傅里叶变换stft

前提：

傅里叶变换是针对平稳信号的，但是很多实际应用中的信号都是非平稳的，如果要计算其傅里叶变换，需要假设其周期无限长，然后对这个无限长的信号做变换分析。但是这种无限长信号分析在实际中并不能实现，我们只能对有限长的局部信号进行分析。

解决：

1、信号在这段时间内是平稳的；

2、在局部平稳信号中截取包含至少1个完整的周期。

对于非周期函数，可以看成是由某个周期函数转换而来的。截取非平稳信号的局部平稳信号，然后采用傅里叶变换进行频域分析，这种分析方法称为短时傅里叶变化。         

存在潜在问题：

如下图所示，若周期截断，则FFT频谱为单一谱线。若为非周期截断，则频谱出现拖尾，如图中部所示，可以看出泄漏很严重。

周期截断，对频率没啥影响，但是非周期截断频率发生很大变化，非周期截断，然后加窗函数，要比原始函数的频率稳定的多。（图片来自网上）

python信号分析实战（stft）：

官网的api：scipy官网api

scipy.signal.stft函数介绍：

scipy.signal.stft(x,fs=1.0,window='hann',nperseg=256,noverlap=None,nfft=None,
detrend=False,return_onesided=True,boundary='zeros',padded=True,axis=-1)

参数说明：

x: 数组类型

时间序列测量值

fs: 浮点型，可选

x时间序列的采样频率，默认1.0。

window ： 字符串或元祖或数组，可选

​需要使用的窗。如果window是一个字符串或元组，则传递给它get_window以生成窗口值，默认情况下使用DFT。有关get_window窗口和所需参数的列表，请参阅。如果window是数组类型，直接以其为窗，其长度必须是nperseg。默认为Hann窗。
nperseg ： int，可选

每个段的长度。默认为256。

noverlap ： int，可选

段之间重叠的点数。如果没有，noverlap = nperseg // 2 .默认为无。如有提前指定时，必须满足COLA约束。

nfft ： int，可选

如果需要零填充FFT，则为使用FFT的长度。如果为 None，则FFT长度为nperseg。默认为无。

detrend ： str或function或False，可选

指定如何去除每个段的趋势。如果detrend是字符串，则将其作为类型参数传递给detrend 函数。如果它是一个函数，它需要一个段并返回一个去趋势段。如果detrend为False，则不进行去除趋势。默认为False。

return_onesided ： bool，可选

如果为True，则返回实际数据的单侧频谱。如果 False返回双侧频谱。默认为 True。请注意，对于复杂数据，始终返回双侧频谱。

boundary ： str或None，可选

指定输入信号是否在两端扩展，以及如何生成新值，以使第一个窗口段在第一个输入点上居中。这具有当所采用的窗函数从零开始时能够重建第一输入点的益处。有效选项是['even', 'odd', 'constant', 'zeros', None].默认为‘zeros’,对于补零操作[1, 2, 3, 4]变成[0, 1, 2, 3, 4, 0] 当nperseg=3.

填充： bool，可选

指定输入信号在末尾是否填充零以使信号精确地拟合为整数个窗口段，以便所有信号都包含在输出中。默认为True。填充发生在边界扩展之后，如果边界不是None，则填充为True，默认情况下也是如此。

axis ： int，可选

计算STFT的轴; 默认值超过最后一个轴（即axis=-1）。

返回值

f ： ndarray

采样频率。

t ： ndarray

时间。

Zxx ： ndarray

x的 STFT 。默认情况下，Zxx的最后一个轴对应于段时间。

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

# 10e3=10000
fs = 10e3
# 1e5=100000，e后面表示10的乘方数
N = 1e5
# 获取矩阵元素的平方根[[1,4],[4,9]]=>[[1,2],[2,3]]
amp = 2 * np.sqrt(2)
noise_power = 0.01 * fs / 2
# 产生一个N长度大小的array
time = np.arange(N)/float(fs)
# numpy.time打印给定弧度时间的余弦值，这的结果是一个aray
mod = 500 * np.cos(2*np.pi*0.25*time)
# 3e3=3000
# numpy.time打印给定弧度时间的正弦值，这的结果是一个aray
carrier = amp * np.sin(2*np.pi*3e3*time+mod)
# np.random.normal()的意思是一个正态分布，normal这里是正态的意思
# 参数loc(float)：正态分布的均值，对应着这个分布的中心。loc=0说明这一个以Y轴为对称轴的正态分布，
# 参数scale(float)：正态分布的标准差，对应分布的宽度，scale越大，正态分布的曲线越矮胖，scale越小，曲线越高瘦。
# 参数size(int 或者整数元组)：输出的值赋在shape里，默认为None。
noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
# 表示指数e的多少次方
noise *= np.exp(-time/5)
# 产生x的一个噪声
x = carrier + noise

# 计算并绘制STFT的大小
f, t, Zxx = signal.stft(x, fs, nperseg=1000)
plt.pcolormesh(t, f, np.abs(Zxx), vmin = 0, vmax = amp)
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

二、短时傅里叶逆变换 istft

scipy.signal.istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, input_onesided=True, boundary=True, time_axis=- 1, freq_axis=- 2)

参数解释：

Zxx： array_like

要重建的信号的 STFT。如果传递的是纯实数数组，它将被强制转换为复杂数据类型。

fs： 浮点数，可选

时间序列的采样频率。默认为 1.0。

window： str 或 tuple 或 数组，可选

想要使用的窗口。如果窗户是一个字符串或元组，它被传递给scipy.signal.get_window生成窗口值，默认为DFT-even。看scipy.signal.get_window获取窗口列表和所需参数。如果窗户是数组，它将直接用作窗口，其长度必须为nperseg。默认为 Hann 窗口。必须与用于生成 STFT 的窗口匹配以实现忠实的反演。

nperseg： int 可选

每个 STFT 段对应的数据点数。如果每段的数据点数是奇数，或者如果 STFT 是通过填充的，则必须指定此参数nfft > nperseg.如果None，值取决于形状Zxx和input_onesided.如果input_onesided是True,nperseg=2*(Zxx.shape[freq_axis] - 1).否则，nperseg=Zxx.shape[freq_axis].默认为None.

noverlap： int 可选

段之间重叠的点数。如果没有，则为段长度的一半。默认为无。指定时，必须满足 COLA 约束(请参阅下面的注释)，并且应与用于生成 STFT 的参数匹配。默认为无。

nfft： int 可选

每个 STFT 段对应的 FFT 点数。如果 STFT 通过以下方式填充，则必须指定此参数nfft > nperseg.如果None, 默认值与nperseg，上面详述，但有一个异常：如果input_onesided是真的并且nperseg==2*Zxx.shape[freq_axis] - 1,nfft也具有该值。这种情况允许使用 odd-length 未填充的 STFT 正确反转nfft=None.默认为None.

input_onesided： 布尔型，可选

如果True，将输入数组解释为one-sided FFT，例如由scipy.signal.stft和return_onesided=True和numpy.fft.rfft.如果False，将输入解释为 two-sided FFT。默认为True.

boundary： 布尔型，可选

指定输入信号是否在其边界处通过提供非扩展None boundary参数scipy.signal.stft.默认为True.

time_axis： int 可选

STFT的时间段所在的位置；默认值为最后一个轴(即 axis=-1 )。

freq_axis： int 可选

STFT的频率轴所在的位置；默认值为倒数第二个轴(即 axis=-2 )。

返回：

t： ndarray

输出数据时间数组。

x： ndarray

Zxx 的 iSTFT。

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

rng = np.random.default_rng()

# 生成一个测试信号，一个 50Hz 的 2 Vrms 正弦波，被 1024 Hz 采样的 0.001 V**2/Hz 白噪声破坏。
# 采样频率1024
fs = 1024
# 数据长度10*1024
N = 10*fs
#窗口大小512
nperseg = 512
# 2*2的开方
amp = 2 * np.sqrt(2)
noise_power = 0.001 * fs / 2
# 采样点对应的时间
time = np.arange(N) / float(fs)
#math.pi == np.pi == scipy.pi   True（输出“3.141592653589793”），求一个弧度值的正弦值，长度和时间点对应
carrier = amp * np.sin(2*np.pi*50*time)
# 随即产生正太分布的噪声，长度和time的时间点一致
noise = rng.normal(scale=np.sqrt(noise_power),
                  size=time.shape)
# 输出的信号就是正弦值+噪声点的混合
x = carrier + noise

# 计算 STFT，并绘制其大小，短时傅里叶变换
f, t, Zxx = signal.stft(x, fs=fs, nperseg=nperseg)
plt.figure()
plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud')
plt.ylim([f[1], f[-1]])
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.yscale('log')
plt.show()

对变换的频率进行滤波，转换小于或者等于10%的幅度 为0幅度。然后进行傅里叶逆变换，逆变换后的波形和原始波形对比。

# 将小于或等于载波幅度 10% 的分量归零，然后通过逆 STFT 转换回时间序列
Zxx = np.where(np.abs(Zxx) >= amp/10, Zxx, 0)
_, xrec = signal.istft(Zxx, fs)

# 将清洁后的信号与原始和真实的载波信号进行比较。
plt.figure()
plt.plot(time, x, time, xrec, time, carrier)
plt.xlim([2, 2.1])
plt.xlabel('Time [sec]')
plt.ylabel('Signal')
plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
plt.show()

 

小知识：

 plt.subplot(221)表示将整个图像窗口分为2行2列, 当前位置为1.

plt.xlim() 显示的是x轴的作图范围

同时plt.ylim() 显示的是y轴的作图范围

而 plt.xticks() 表达的是x轴的刻度内容的范围

缩小x的范围： 

plt.figure()
plt.plot(time, x, time, xrec, time, carrier)
plt.xlim([0, 0.1])
plt.xlabel('Time [sec]')
plt.ylabel('Signal')
plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
plt.show()

总结：利用傅里叶短时变换和逆变换，可以过滤信号中我们不需要的部分，对信号进行改造。