文章目录
一、短时傅里叶变换（Short-Time Fourier Transform，STFT）
1. 基本思想
把信号划分为许多小的时间间隔，用傅里叶变换分析每一个时间间隔，以便确定该时间间隔存在的频率。
2. 定义
给定一个时间宽度很短的窗函数 η ( t ) \eta(t) η(t)，让窗滑动，则信号 s ( t ) s(t) s(t)的短时傅里叶变换（STFT）可以定义为
S T F T s ( t , f ) = ∫ − ∞ ∞ s ( u ) η ∗ ( u − t ) e − j 2 π f u d u (1) STFT_s(t,f)=\int_{-∞}^{∞}s(u)\eta^*(u-t)e^{-j2\pi fu}\,du \tag{1} STFTs(t,f)=∫−∞∞s(u)η∗(u−t)e−j2πfudu(1)
式中， ∗ ^* ∗表示复共轭。
由式可见，正是由于窗函数 η ( t ) \eta(t) η(t)的存在，使得短时傅里叶变换成为时间和频率的二维函数，它将信号的时域和频域联系起来，可以据此对信号进行时频分析。
定义式表明，信号s(u)在时间 t 0 t_0 t0处的短时傅里叶变换就是信号乘上一个短窗函数 η ( u − t 0 ) \eta(u-t_0) η(u−t0)，即取出信号在分析时间点 t 0 t_0 t0 附近的一个切片，所以短时傅里叶变换 S T F T s ( t , f ) STFT_s(t,f) STFTs(t,f)可以理解为信号s(u)在时间 t 附近的傅里叶变换，即“局部频谱”。
3. 短时傅里叶变换的完全重构条件
为了使短时傅里叶变换真正成为一种有实际价值的非平稳信号分析工具，信号s(t)应该能够由 S T F T s ( t , f ) STFT_s(t,f) STFTs(t,f)完全重构出来。
记重构窗函数为g(t)，则要求窗函数满足以下条件：
∫ − ∞ ∞ η ∗ ( t ) g ( t ) d t = 1 (2) \int_{-∞}^{∞}\eta^*(t)g(t)\,dt=1 \tag{2} ∫−∞∞η∗(t)g(t)dt=1(2)
证明如下：假设重构公式为
β ( λ ) = ∫ − ∞ ∞ ∫ − ∞ ∞ S T F T s ( t , f ) g ( λ − t ) e j 2 π f λ d t d f (3) β(\lambda)=\int_{-∞}^{∞}\int_{-∞}^{∞}STFT_s(t,f)g(\lambda-t)e^{j2\pi f \lambda}\,dtdf \tag{3} β(λ)=∫−∞∞∫−∞∞STFTs(t,f)g(λ−t)ej2πfλdtdf(3)
将式(1)代入式(3)，可以得到
β ( λ ) = ∫ − ∞ ∞ ∫ − ∞ ∞ [ ∫ − ∞ ∞ s ( u ) η ∗ ( u − t ) e − j 2 π f u d u ] g ( λ − t ) e j 2 π f λ d t d f = ∫ − ∞ ∞ ∫ − ∞ ∞ [ ∫ − ∞ ∞ e − j 2 π f ( u − λ ) d f ] s ( u ) η ∗ ( u − t ) g ( λ − t ) d u d t = ∫ − ∞ ∞ ∫ − ∞ ∞ s ( u ) η ∗ ( u − t ) g ( λ − t ) δ ( u − λ ) d u d t = s ( λ ) ∫ − ∞ ∞ η ∗ ( λ − t ) g ( λ − t ) d t = s ( λ ) ∫ − ∞ ∞ η ∗ ( t ) g ( t ) d t \begin{aligned} β(\lambda)&=\int_{-∞}^{∞} \int_{-∞}^{∞} \Big [\int_{-∞}^{∞}s(u) \eta^*(u-t)e^{-j2\pi fu}\,du \Big] g(\lambda-t)e^{j2\pi f\lambda}\,dtdf\\ &=\int_{-∞}^{∞} \int_{-∞}^{∞} \Big [\int_{-∞}^{∞} e^{-j2\pi f(u-\lambda)}\,df \Big]s(u)\eta^*(u-t) g(\lambda-t)\,dudt\\ &=\int_{-∞}^{∞} \int_{-∞}^{∞}s(u)\eta^*(u-t) g(\lambda-t)\delta(u-\lambda)\,dudt\\ &=s(\lambda)\int_{-∞}^{∞} \eta^*(\lambda-t) g(\lambda-t)\,dt\\ &=s(\lambda)\int_{-∞}^{∞} \eta^*(t) g(t)\,dt\\ \end{aligned} β(λ)=∫−∞∞∫−∞∞[∫−∞∞s(u)η∗(u−t)e−j2πfudu]g(λ−t)ej2πfλdtdf=∫−∞∞∫−∞∞[∫−∞∞e−j2πf(u−λ)df]s(u)η∗(u−t)g(λ−t)dudt=∫−∞∞∫−∞∞s(u)η∗(u−t)g(λ−t)δ(u−λ)dudt=s(λ)∫−∞∞η∗(λ−t)g(λ−t)dt=s(λ)∫−∞∞η∗(t)g(t)dt
因此，为了完全重构信号，即使得 β ( λ ) = s ( λ ) β(\lambda)=s(\lambda) β(λ)=s(λ)，则必须满足条件 ∫ − ∞ ∞ η ∗ ( t ) g ( t ) d t = 1 \int_{-∞}^{∞}\eta^*(t)g(t)\,dt=1 ∫−∞∞η∗(t)g(t)dt=1
证毕。
完全重构条件是一个很宽松的约束条件，对于一给定的分析窗函数 η ( t ) \eta(t) η(t)，满足条件式(2)的窗函数 g ( t ) g(t) g(t)有无穷种可能的选择。
特别地，取 g ( t ) = η ( t ) g(t)=\eta(t) g(t)=η(t)，此时对应的完全重构条件式(2)变为
∫ − ∞ ∞ η ∗ ( t ) g ( t ) d t = ∫ − ∞ ∞ ∣ η ( t ) ∣ 2 d t = 1 (4) \int_{-∞}^{∞}\eta^*(t)g(t)\,dt=\int_{-∞}^{∞}|\eta(t)|^2\,dt=1 \tag{4} ∫−∞∞η∗(t)g(t)dt=∫−∞∞∣η(t)∣2dt=1(4)
即所谓能量归一化，此时式(3)可写成
s ( λ ) = ∫ − ∞ ∞ ∫ − ∞ ∞ S T F T s ( t , f ) η ( λ − t ) e j 2 π f λ d t d f (5) s(\lambda)=\int_{-∞}^{∞}\int_{-∞}^{∞}STFT_s(t,f)\eta(\lambda-t)e^{j2\pi f \lambda}\,dtdf \tag{5} s(λ)=∫−∞∞∫−∞∞STFTs(t,f)η(λ−t)ej2πfλdtdf(5)
式(5)称之为广义短时傅里叶反变换。
与正、反傅里叶变换不同的是，短时傅里叶变换是一维变换，而广义短时傅里叶反变换是二维变换。
4. 离散形式
在实际应用中，需要对 S T F T s ( t , f ) STFT_s(t,f) STFTs(t,f)进行离散化处理，为此在时频面上等间隔时频网络点 ( m ⋅ Δ t ， n ⋅ Δ f ) (m\cdot\Delta t，n\cdot\Delta f) (m⋅Δt，n⋅Δf)处采样，其中 Δ t \Delta t Δt， Δ f \Delta f Δf分别表示时间变量和频率变量的采样间隔，令s(k)表示信号s(t)的离散形式，则短时傅里叶变换的离散形式为
S T F T s ( m , n ) = ∑ k = − ∞ ∞ z ( k ) η ( k Δ t − m Δ t ) e − j 2 π ( n Δ f ) k (6) STFT_s(m,n)=\sum_{k=-∞}^{∞}z(k)\eta(k \Delta t-m\Delta t)e^{-j2\pi(n\Delta f)k} \tag{6} STFTs(m,n)=k=−∞∑∞z(k)η(kΔt−mΔt)e−j2π(nΔf)k(6)
同样，式(5)可以离散为
s ( k ) = ∑ m = − ∞ ∞ ∑ n = − ∞ ∞ S T F T s ( m , n ) η ( k Δ t − m Δ t ) e j 2 π ( n Δ f ) k (7) s(k)=\sum_{m=-∞}^{∞}\sum_{n=-∞}^{∞}STFT_s(m,n)\eta(k \Delta t-m\Delta t)e^{j2\pi(n\Delta f)k} \tag{7} s(k)=m=−∞∑∞n=−∞∑∞STFTs(m,n)η(kΔt−mΔt)ej2π(nΔf)k(7)
二、MATLAB实现
1. 编写程序绘制时频图
clear all;
close all;
clc;
%%%%%%%%%% 应用短时傅里叶变换STFT原理绘制信号时频图 %%%%%%%%%
%信号参数
TimeWidth=10e-6; % 信号时宽
BandWidth=30e6; % 信号带宽
Fs=2*BandWidth; % 采样频率
F0=3/4*Fs; % 中频频率
%信号定义
N=fix(TimeWidth*Fs); % 信号长度，fix函数向零方向取整
t=(-N/2:N/2-1)/Fs;
st=cos(2*pi*F0*t+pi*(BandWidth/TimeWidth)*t.^2); % 线性调频信号余弦表达式
%STFT参数设置
win=hamming(128); % 窗的类型
wlen=length(win); % 窗长，即做一次傅里叶变换的长度
noverlap=120; % 重叠的长度
step=wlen-noverlap; % 每次移动的步进
step_nums=fix((N-wlen)/step)+1; % 需要移动的步数，进行短时傅里叶变换的个数
%拆分信号
for i=1:step_nums
ts=(1:wlen)+(i-1)*step;
x(:,i)=st(ts).*win'; % 提取第i段信号数据,加窗，存储在x矩阵的第i列
S(:,i)=fftshift(abs(fft(x(:,i)))); % 第i段信号的傅里叶变换，存储在S矩阵的第i列
end
tt= (wlen/2:step:wlen/2+(step_nums-1)*step)/Fs; % 取每个窗口中点的时间点
f = (-N/2:N/2-1) * (Fs/N); % 频率范围
%绘制时频图
imagesc(tt*1e6,f*1e-6,S); % 绘制图像，并使图像位于指定的区域
xlabel('time/us'); ylabel('frequency/MHz'); title('线性调频信号的时频图');
2. 调用函数绘制时频图
2.1 spectrogram
2.2 stft
clear all;
close all;
clc;
%%%%%%%%%% 调用MATLAB自带函数spectrogram\stft绘制信号时频图 %%%%%%%%%%
% 信号参数
TimeWidth=10e-6; %信号时宽
BandWidth=30e6; %信号带宽
Fs=60e6; %采样频率
F0=3/4*Fs; %中频频率
% 信号定义
N=fix(TimeWidth*Fs); %采样点数，fix函数向零方向取整
t=(-N/2:N/2-1)/Fs;
st=cos(2*pi*F0*t+pi*(BandWidth/TimeWidth)*t.^2); %线性调频信号余弦表达式
%%%%%%%%%% 1.调用函数spectrogram绘制信号时频图 %%%%%%%%%%
% 函数参数设置
win=hamming(128);
noverlap=120;
nfft=2^nextpow2(length(win));
figure;
[s1,f1,t1]=spectrogram(st,win,noverlap,nfft,Fs,'yaxis'); %注意：spectrogram输出的是单边谱
imagesc(t1*1e6,f1*1e-6,abs(s1));
xlabel('time/us'); ylabel('frequency/MHz'); title('spectrogram');
%%%%%%%%%% 2.调用函数stft绘制信号时频图 %%%%%%%%%%
[s2,f2,t2]=stft(st,Fs,'Window',hamming(128),'OverlapLength',120);
figure;
imagesc(t2*1e6,f2*1e-6,abs(s2));
xlabel('time/us'); ylabel('frequency/MHz'); title('stft');
三、窗口长度与分辨率
短时傅里叶变换的分辨率和窗口长度n的选取，平移的步长step都有关系。
当窗口大到整个信号长度时，短时傅里叶变换就退化为傅里叶变换，没有时间维度。当窗口小到单个采样点时，分析也就退化为时域分析，无法做频域分析。
注意 S 的尺寸随 wlen 的变换，会发现：
- 一个好的时间分辨率需要一个短的窗函数，然而，窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差；
- 一个高的频率分辨率需要一个长的窗函数，然而，窗太宽，时域上又不够精细，时间分辨率低。
根据Heisenberg不确定性原理，时间分辨率与频率分辨率不能同时任意小，它们的乘积受到一定值的限制。要提高时间分辨率就要降低频率分辨率，反之亦然。
四、短时傅里叶变换的局限性
短时傅里叶变换虽然可以描述某一局部时间段上的频率信息，但是其时、频域的分辨率 Δ t ， Δ f \Delta t，\Delta f Δt，Δf不随时间 t 和频率 f 的变化而变化。
实际应用中一旦窗函数选定，则时间分辨率与频率分辨率确定。对于我们要分析的非平稳信号而言，也许某一小时段上是以高频信息为主，我们希望用短时间窗进行分析；而在某一长时间段上是一些低频信息，我们希望用一个长时间窗进行分析。
因此具有不变窗的短时傅里叶变换，更适合应用在准稳态信号分析的场合。对于一个时变的非平稳信号，利用短时傅里叶变换方法难以找到一个合适的时间窗口来适应于不同的时间段。
扫一扫
7万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
