传统的信号处理主要基于信号时域和频域来进行，信号的时域信息表征了信号幅度随时间的变化关系，信号的频域信息则揭示了信号频率和能量之间的关系。电磁信号通常都是典型的非平稳信号，只进行时域角度的分析则无法获取频域的信息，而只分析了全局频谱的傅里叶变换不仅无法对非平稳信号进行准确有效的描述，还丢失了信号的时域信息。从 20 世纪 40 年代起，科研人员就开始将数学领域的理论应用到信号处理中，运用时频分析方法对电磁信号进行分析和处理，可以同时表征信号的时域和频域信息，时频分析已成为现代数字信号处理中一种常见的信号特征提取方法。
一、短时傅里叶变换
1.1 理论
对于平稳信号而言，傅里叶变换（Fourier Transform, FT）作为一种分析和处理的重要方法，被广泛应用于数学、信号处理、热力学分析和工程控制论等领域。但对于自然界中广泛存在的非平稳信号，傅里叶变换只能得到其全局的频谱信息而丢失在局部时间内的频率变化信息，无法准确地表征非平稳信号的特性。连续时间信号s(t)s(t)s(t)的傅里叶变换定义如下：
S(f)=F[s(t)]=∫−∞∞s(t)e−j2πftdt
S(f)=\mathcal{F}[s(t)]=\int_{-\infty}^{\infty}s(t)e^{-j2\pi ft}\mathrm{d}t
S(f)=F[s(t)]=∫−∞∞s(t)e−j2πftdt
1946 年 Gabor 为了解决由于傅里叶变换对信号频谱表征的全局性而无法体现出电磁信号频率随时间变化规律的问题，提出了短时傅里叶变换(Short Time Fourier Transform, STFT)的概念，其基本思想是在信号傅里叶变换前乘上一个固定长度的矩形窗函数 ，从而实现信号在局部时间里的伪平稳，然后通过矩形窗在时间轴上的移动对信号逐段进行傅里叶变换，从而得到信号的时变特性。
连续时间信号s(t)s(t)s(t)的 STFT 定义如下：
STFT(t,f)=∫−∞∞s(τ)h(τ−t)e−j2πfτdτ
STFT(t,f)=\int_{-\infty}^{\infty}s(\tau)h(\tau-t)e^{-j2\pi f\tau}\mathrm{d}\tau
STFT(t,f)=∫−∞∞s(τ)h(τ−t)e−j2πfτdτ
离散时间信号s[n]s[n]s[n]的 STFT 定义如下：
STFT[m,k]=∑ns[n]ω[n−k]e−j2πNmnm=0,1,...,N−1
STFT[m,k]=\sum_ns[n]\omega[n-k]e^{-j\frac{2\pi}{N}mn} \quad\quad\quad m=0,1,...,N-1
STFT[m,k]=n∑s[n]ω[n−k]e−jN2πmnm=0,1,...,N−1
时间窗函数的宽度决定了信号 STFT 的时间分辨率，窗函数宽度TpT_pTp越小，时间分辨率越高，二者关系表达式为：
Tp=NT=Nfsample
T_{p}=NT=\frac{N}{f_{sample}}
Tp=NT=fsampleN
频率分辨率是度量信号 STFT 效果的重要参数，是指分辨信号频谱中最小相邻谱峰间隔的能力，Δfc\Delta f_cΔfc越小频率分辨率越高，与时间分辨率的关系如下：
Δfc=1Tp
\Delta f_c=\frac1{T_p}
Δfc=Tp1
短时傅里叶变换计算方法相对简单，是一种联合的信号线性时频分析方法。利用固定长度的矩形窗函数对信号进行时频分析，可以避免高次型非平稳信号分析方法中出现的交叉项干扰问题，适用于信号的多分量分析。
1.2 MATLAB实现
在MATLAB中，可以使用stft函数来实现STFT。这个函数可以自动处理窗函数、窗口长度、步长等参数，也根据自定义参数来适应不同的信号分析需求。
clc
clear
close all
set(0,'defaultfigurecolor','w') %设置图窗背景为白色
% 生成一个非平稳信号
fs = 1000; % 采样频率 (Hz)
t = 0:1/fs:2; % 时间向量
x = cos(2*pi*50*t) + cos(2*pi*120*t).*(t > 0.5 & t < 1.5); % 混合信号
% 设置STFT参数
window = hamming(128); % 窗函数 (128点汉明窗)
noverlap = 64; % 窗口重叠样本数
nfft = 256; % FFT点数
% 计算STFT
[stft_result, f, t_stft] = stft(x, fs, 'Window', window, 'OverlapLength', noverlap, 'FFTLength', nfft);
% 绘制STFT的时频图
figure;
surf(t_stft, f, abs(stft_result), 'EdgeColor', 'none');
axis tight;
view(0, 90);
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('STFT Magnitude');
colorbar;
参数说明：
fs：采样频率。window：窗函数，用于选择信号的窗口长度和形状。noverlap：窗口重叠样本数，通常取窗口长度的一半。nfft：FFT的点数，决定频率分辨率。stft_result：STFT计算的结果，是一个复数矩阵，取abs(stft_result)可以得到幅度谱。f 和 t_stft：分别对应频率和时间轴的坐标。
生成的时频图像如图所示
二、小波变换
2.1 理论
通过在时间上将非平稳信号分割成许多小的时间间隔而实现非平稳信号的局部伪平稳，然后在对每个伪平稳局部信号进行傅里叶变换的基础上乘以矩形窗函数得到信号的短时傅里叶变换。虽然短时傅里叶变换克服了傅里叶变换的频谱分布全局性的特点，但是由于其选取的矩形窗函数的宽度是固定的，无法同时兼顾较高的时间分辨率和频率分辨率，较长时间宽度通过牺牲时间分辨率得到理想的频率分辨率；较短时间宽度则通过牺牲频率分辨率来得到理想的时间分辨率。
连续小波变换（Continue Wavelet Transform, CWT）继承和发展了短时傅里叶对非平稳信号处理采取的局部分割的思想，同时克服短时傅里叶窗函数固定宽度的缺点，为信号进行时频变换提供了一个随频率改变的时间窗口。信号连续小波变换的定义如下：
CWT(x,y)=1x∫−∞∞s(t)φ(t−yb)dt
CWT(x,y)=\frac1{\sqrt{x}}\int_{-\infty}^\infty s\left(t\right)\varphi{\left(\frac{t-y}b\right)}dt
CWT(x,y)=x1∫−∞∞s(t)φ(bt−y)dt
式中的 xxx 和 yyy 分别表示尺度和平移量。尺度 xxx 表征小波基函数的伸缩量，与频率成反比，可以在时频图中表征频率；平移量 yyy 表示小波函数的平移量，对应于时频图中时间的长短。
小波分析在信号处理、图像识别、语音合成和控制论等领域都取得了广泛的应用和巨大的成果。在通信信号处理中，小波变换主要用于分离噪声、提取弱信号、时频分析、滤波和信号识别等。
2.2 MATLAB实现
在MATLAB中可以使用 cwt 函数来计算并可视化小波变换。
% 生成一个非平稳信号
fs = 1000; % 采样频率 (Hz)
t = 0:1/fs:2; % 时间向量
x = cos(2*pi*50*t) + cos(2*pi*120*t).*(t > 0.5 & t < 1.5); % 混合信号
% 使用小波变换计算信号的时频表示
[cwt_result, f] = cwt(x, fs);
% 绘制小波变换的时频图
figure;
surface(t, f, abs(cwt_result), 'EdgeColor', 'none');
axis tight;
view(0, 90);
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('CWT Magnitude');
colorbar;
参数说明：
fs：采样频率，用于小波变换的采样频率参数。cwt_result：小波变换结果，得到的是一个复数矩阵，取abs(cwt_result)可以获得幅度谱。f：小波变换的频率轴，与信号的频率分量对应。
生成的时频图像如图所示
三、魏格纳-维利分布
3.1 理论
1932年美国物理学家Eugene Wigner提出了Wigner-Ville分布（WVD）的概念，Vile于1948年将WVD应用到信号处理领域中。连续时间信号 的Wigner-Ville分布可以通过对信号自身的瞬时自相关函数进行傅里叶变换得到，信号的瞬时自相关表达式如下：
rs(t,τ)=s(t+τ2)s∗(t−τ2)
r_s(t,\tau)=s\Big(t+\frac{\tau}{2}\Big)s^*\Big(t-\frac{\tau}{2}\Big)
rs(t,τ)=s(t+2τ)s∗(t−2τ)
其WVD的定义式如下：
WVD(t,f)=∫−∞∞s(t+τ2)s∗(t−τ2)e−j2πfτdτ
WVD\left(t,f\right)=\int_{-\infty}^{\infty}s{\left(t+\frac{\tau}{2}\right)}s^{*}\left(t-\frac{\tau}{2}\right)e^{-j2\pi f\tau}\mathrm{d}\tau
WVD(t,f)=∫−∞∞s(t+2τ)s∗(t−2τ)e−j2πfτdτ
信号的魏格纳-维利分布具有较高的时间分辨率和频率分辨率，并且有良好的时频聚集性，克服了 STFT 的缺点，适合分析非稳态的随机信号。
3.2 MATLAB实现
在MATLAB中，魏格纳-维利分布（Wigner-Ville Distribution，WVD）没有内置的直接函数，但可以使用信号处理工具箱中的一些方法来实现。
% 生成一个非平稳信号
fs = 1000; % 采样频率 (Hz)
t = 0:1/fs:2; % 时间向量
x = cos(2*pi*50*t) + cos(2*pi*120*t).*(t > 0.5 & t < 1.5); % 混合信号
% 魏格纳-维利分布计算
n = length(x);
wvd_result = zeros(n, n); % 初始化WVD矩阵
for tau = -n/2+1:n/2 % 使用整数偏移量
% 计算时延信号和原信号的乘积
shifted_x = circshift(x, round(tau)); % 确保偏移量为整数
wvd = x .* conj(shifted_x); % 时域自相关函数
% 傅里叶变换得到时频分布
wvd_result(:, mod(tau + n/2 - 1, n) + 1) = fft(wvd);
end
% 只保留前半部分频率
f = linspace(0, fs/2, n/2);
wvd_result = abs(wvd_result(1:n/2, :));
% 绘制WVD的时频图
figure;
imagesc(t, f, wvd_result);
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('Wigner-Ville Distribution');
colorbar;
参数说明
fs：采样频率。t：时间轴，用于绘制结果的横坐标。wvd_result：魏格纳-维利分布计算的结果矩阵，取其绝对值以获得幅度谱。f：频率轴，对应的是时频分布的频率分量。
生成的时频图像如图所示
四、希尔伯特-黄变换
4.1 理论
1998 年，Norden E. Huang 等人提出了经验模态分解（Empirical Mode Decomposition, EMD）方法，并引入了希尔伯特谱的概念和分析方法，这一方法后来被命名为 Hilbert-Huang Transform，简称 HHT，即希尔伯特-黄变换。HHT 处理非平稳信号的基本思路是：首先利用经验模态分解方法将给定的信号分解为若干固有模态函数（Intrinsic Mode Function, IMF），然后对每一个固有模态函数进行希尔伯特变换得到相应的希尔伯特谱，最后汇总所有固有模态函数的希尔伯特谱就得到原始信号的希尔伯特谱。
希尔伯特-黄变换主要包括两个过程，经验模态分解和希尔伯特变换，经验模态分解方法过程如下：
（1）分别提取待分析信号s(t)s(t)s(t)的所有局部极大值和极小值点，再分别拟合极大值点、极小值点形成极大值和极小值包络线，得到极值包络线均值。
m1=emax(t)+emin(t)2
m_1=\frac{e_{max}(t)+e_{min}(t)}2
m1=2emax(t)+emin(t)
（2）信号s(t)s(t)s(t)减去极值包络线的均值m1m_1m1得到新信号y1(t)y_1(t)y1(t)
y1(t)=s(t)−m1
y_1(t)=s(t)-m_1
y1(t)=s(t)−m1
得到的新信号y1(t)y_1(t)y1(t)就是固有模态函数的分量形式，然后判断该分量形式是否满足条件。如果y1(t)y_1(t)y1(t)不满足条件，则将y1(t)y_1(t)y1(t)作为原始数据重复步骤（1）和（2），不断筛选一直到y1(t)y_1(t)y1(t)满足固有模态函数的分量条件，记满足条件的y1(t)=c1(t)y_1(t)=c_1(t)y1(t)=c1(t)，则c1(t)c_1(t)c1(t)是信号s(t)s(t)s(t)的第一个固有模态函数分量。
（3）将剩余信号r1(t)=s(t)−c1(t)r_{1}(t)=s(t)-c_{1}(t)r1(t)=s(t)−c1(t)作为原始信号，同上面原理分解 jjj次得到 jjj 阶固有模态函数，余量为rn(t)r_n(t)rn(t)。如果得到的余量rn(t)r_n(t)rn(t)是一个单调信号或mnm_nmn的值足够小时，经验模态分解过程结束。综合上述结果，将各阶固有模态函数及余项rn(t)r_n(t)rn(t)重构获得如下的原始信号：
s(t)=∑i=1nci(t)−cn(t)
s(t)=\sum_{i=1}^nc_i(t)-c_n(t)
s(t)=i=1∑nci(t)−cn(t)
对原信号进行经验模态分解后，再对分解得到的固有模态函数分量进行希尔伯特变换就可以得到每个固有模态函数的瞬时频率，然后将所有分量的瞬时频谱综合即可获得原始信号的希尔伯特谱。
各阶固有模态函数分量进行希尔伯特变换：
H[c(t)]=1πPV∫−∞∞s(τ)t−τdτ
H[c(t)]=\frac{1}{\pi}PV\int_{-\infty}^{\infty}\frac{s(\tau)}{t-\tau}\mathrm{d}\tau
H[c(t)]=π1PV∫−∞∞t−τs(τ)dτ
式中 PVPVPV 表示柯西主值，根据瞬时频率公式的定义：
ω(t)=d[arctan(s^(t)s(t))]dt
\omega(t)=\frac{d\Big[\arctan\left(\frac{\hat{s}(t)}{s(t)}\right)\Big]}{dt}
ω(t)=dtd[arctan(s(t)s^(t))]
可以将原始信号表示为：
s(t)=Re(∑i=1nai(t)ejϕi(t))=Re(∑i=1nai(t)ej∫ωi(t)dt)
s(t)=\text{Re}\left(\sum_{i=1}^{n}a_{i}\left(t\right)e^{j\phi_{i}\left(t\right)}\right)=\text{Re}\left(\sum_{i=1}^{n}a_{i}\left(t\right)e^{j\int\omega_{i}\left(t\right)dt}\right)
s(t)=Re(i=1∑nai(t)ejϕi(t))=Re(i=1∑nai(t)ej∫ωi(t)dt)
展开后就是希尔伯特谱，可以写成：
H(ω,t)=Re(∑i=1nai(t)ej∫ωi(t)dt)
H\left(\omega,t\right)=\mathrm{Re}\left(\sum_{i=1}^{n}a_{i}\left(t\right)e^{j\int\omega_{i}\left(t\right)dt}\right)
H(ω,t)=Re(i=1∑nai(t)ej∫ωi(t)dt)
希尔伯特谱是一种完整的时间、频率和能量的联合谱，其中可以同时看到不同时间段的频率变化情况和能量随着时间与频率的变化而变化的情况。
4.2 MATLAB实现
在MATLAB中实现希尔伯特-黄变换（HHT）可以通过两步完成：首先使用经验模态分解（EMD）将信号分解为若干固有模态函数（IMF），然后对每个IMF应用希尔伯特变换以获得瞬时频率。MATLAB的信号处理工具箱提供了 emd 函数来实现经验模态分解，hilbert 函数则可以用于计算希尔伯特变换。
% 生成一个非平稳信号
fs = 1000; % 采样频率 (Hz)
t = 0:1/fs:2; % 时间向量
x = cos(2*pi*50*t) + cos(2*pi*120*t).*(t > 0.5 & t < 1.5); % 混合信号
% 1. 经验模态分解 (EMD)
imfs = emd(x); % 分解得到IMFs
% 初始化变量以存储每个IMF的瞬时频率
num_imfs = size(imfs, 2); % IMF的数量
inst_freqs = zeros(num_imfs, length(t) - 1); % 用于存储每个IMF的瞬时频率
% 2. 对每个IMF应用希尔伯特变换
for i = 1:num_imfs
% 计算IMF的希尔伯特变换
imf_hilbert = hilbert(imfs(:, i));
% 瞬时频率计算
inst_freqs(i, :) = fs / (2 * pi) * diff(unwrap(angle(imf_hilbert))); % 瞬时频率
end
% 3. 绘制时频图
figure;
for i = 1:num_imfs
subplot(num_imfs, 1, i);
plot(t(2:end), inst_freqs(i, :), 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title(['IMF ' num2str(i) ' Instantaneous Frequency']);
end
sgtitle('Hilbert-Huang Transform (HHT)');
参数说明：
fs：采样频率。imfs：使用EMD分解得到的固有模态函数（IMFs）。inst_freqs：存储每个IMF的瞬时频率。hilbert：对每个IMF应用希尔伯特变换以获取瞬时相位，从而计算瞬时频率。
代码解释：
EMD分解：将信号分解成若干IMFs，适合捕捉信号的本征振荡模式。希尔伯特变换：对每个IMF进行希尔伯特变换，计算出瞬时相位，得到瞬时频率。绘图：逐个绘制IMF的瞬时频率曲线，展示信号的时频特性。
信号的HHT变换过程图像如下