MATLAB 实现相关分析和谱估计

对于没有系统学过《信号与原理》的自己来说，即使已经在用各种信号处理的方法，仍然在每次接触相关的概念和方法时觉得陌生和模糊。本篇文章整理自己学习过程中的理解和备忘，希望对同样的信号分析小白有所帮助。感谢网上各种资源和分享对我学习过程的帮助。

PS 持续更新中……请多多包涵꒦ິ^꒦ິ

一. 基本概念

1.1 什么是信号？

信号–$f(t)$：信号表示，时域，频域（与傅里叶变换），信号调制（与边频带和倒频谱）

信号在数学上即为一个或多个变量的函数，相应的，对函数的处理可以适用于对信号的处理中。
函数的平移、对称和缩放，对应了信号在时域上的基本变换（时移、反转–相关性和卷积会用到，尺度变换–小波变换会用到）。
从函数的周期性出发，我们便有了频域角度来描述信号的初级理解。
（有限个不同周期性的函数叠加，可以在时域得到一个连续的周期函数。通过傅里叶变换来实现时域和频域的相互转换，频域非周期离散，时域周期连续。而时域非周期连续，频域也会是非周期连续的。——）

能量–$E$：$$E=\underset{T\to \infty }{\lim }{\int }_{-T}^T|f(t){|}^{2}dt$$
功率–$P$：$$P=\underset{T\to \infty }{\lim }\frac{1}{2T}{\int }_{-T}^T|f(t){|}^{2}dt$$
功率和能量的概念其实是类比物理中的定义给出的。把信号值视为电压信号，电阻视为 1，电压的平方就对应能量，能量除以时间就对应功率。

确定信号：确定信号分为能量信号和功率信号。
随机信号：随机信号只有功率信号。

功率（有限）信号：如各种周期信号、常值信号、阶跃信号等。无限长的白噪音也是功率信号。功率有限但能量无限。
能量（有限）信号：如各类瞬变信号。所有有限数量的脉冲信号都是能量信号。有限能量，零功率（因为相当于周期无穷大）。

参考资料：
信号&系统 | 信号的概念与常用信号
随机过程（随机信号）
能量信号和功率信号的分别
傅里叶分析之掐死教程（完整版）更新于2014.06.06
傅里叶变换的基本性质

从时域到频域再到复数域！
傅里叶变换则可以让微分和积分在频域中变为乘法和除法

1.2 信号频域分析方法的理解

一般来说，频谱分析指的是将信号做傅里叶变换从而进行分析。（“谱”意味着自变量为频率，也就是在频域进行分析。）
完整的频域分析应该得到被测信号包含的频率成分，还有每个频率成分的幅度和相位关系。即信号功率谱和相位谱的分析。

FFT（快速傅里叶变换，离散傅里叶变换的快速算法）的假设条件：根据狄里赫利条件，只有能量信号才可以做 FFT 得到频谱。
对于能量信号，能量谱（密度）是原信号傅里叶变换的平方。

功率谱（密度）是针对随机信号分析提出的概念。（随机信号一定是功率信号）
对于功率信号，先求自相关，再做傅里叶变换，物理意义上就是功率谱。对应间接法求功率谱（见1.4）。
功率信号也可以先进行有限截断（变为能量信号），这样才可以做 FFT。再对 FFT 幅值做平方，对应直接法求功率谱（见1.4）。

总结一下¹：

谱密度：自谱和互谱
谱分析：互谱密度
谱密度估计：为什么要估计。非平稳信号来着？

参考资料：
信号频域分析方法的理解（频谱、能量谱、功率谱、倒频谱、小波分析）
如何通俗地讲解傅立叶分析和小波分析间的关系？

1.3 信号的（时滞）相关性

（不同类别的信号的自相关函数）
（对相关函数的解释，比如周期性、区分出不同类别的信号、局域相关）

相关分析：时域卷积，频域乘积
幅值调制：时域乘积，频域卷积

相关函数:

相关系数:

互相关只看相位信息，与振幅无关。
参考资料中有简单的代码示例，可以发现，只改变信号的幅值大小得到的相关系数结果不会有变化。

应用：

1. 时移测距–雷达；时移→相移 $fdt$
2. 噪声信号的自相关函数会随着时移的增大而迅速减小
   这样的特点，在做周期信号测量的时候，可以用自相关或互相关将随机噪声滤除掉，在工程上叫做相关滤波

如果有非白噪声的其他特征频率，先提取目标频率再进行相关分析是有必要的吗？

参考资料：
利用matlab函数xcorr对信号进行相关分析（包含为何加上无偏估计参数）
信号的相关分析 | 相关系数+相关函数+相关定理
自相关与互相关
用matlab计算波形互相关系数(Cross-Correlation)

1.4 功率谱密度 PSD 估计方法

功率谱密度 PSD：单位频带内的信号功率。与 FFT 得到的频谱相比，功率谱丢失了频谱的相位信息。
谱密度估计实现方法：直接法（周期图法）、间接法（自相关函数法）；对序列或估计的自相关函数加窗处理，前者称为数据窗，后者称为滞后窗。[2]

• 直接法：随机序列 $x(n)$ 的 $N$ 个观测数据视为一能量有限的序列，直接计算 $x(n)$ 的离散傅立叶变换，得 $X(k)$ ，然后再取其幅值的平方，并除以 $N$ ，作为序列 $x(n)$ 真实功率谱的估计。
• 间接法：先由序列 $x(n)$ 估计出自相关函数 $R(n)$ ，然后对 $R(n)$ 进行傅立叶变换，便得到 $x(n)$ 的功率谱估计。功率谱与自相关函数是一个傅氏变换对。-- 维纳-辛钦定理

从原理上讲似乎没什么区别，从MATLAB仿真结果上来看，相关函数法对噪声的抑制效果更好，图线更平滑。²
（谱估计：图示和公式证明。）

方法改进：
Bartlett 平均周期图的方法是将 $N$ 点的有限长序列 $x(n)$ 分段求周期图再平均。
Welch 法对 Bartlett 法进行了两方面的修正：Welch’s Overlapped Segment Averaging Spectral Estimation

• 一是选择适当的窗函数 $w(n)$ ，并在周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负（？）；
• 二是在分段时，可使各段之间有重叠，这样会使方差减小[3]。

psepctrum/Signal Analyzer APP 中，如果通过分析整个段得出的分辨率无法实现，APP 将计算 Welch 周期图。[6]

可参考资料：
频域特征值提取的MATLAB代码实现（频谱、功率谱、倒频谱）

1.5 窗函数和频谱泄露

（窗函数：图示和公式表示。）

（如何选择窗函数的形状和长度？）
（加窗和信号有限长的关系？）
频率分辨率和谱泄露分别对应加窗的频率变换中的主瓣和旁瓣宽度，分辨率越好，谱泄露越高。比如，矩形窗有最窄主带和最高的边带(?)，可以分辨两个距离很近的且能量近似的频率成分，但无法检测出较弱的信号。[5]

当 β=5.7 | kaiser(l,beta) 即 l=1-(β/40)=0.8575 | pspetrum(…,‘Leakage’,l) 时，Kaiser 窗和 Hann 窗类似；当 l=1 即 β=0 Kaiser 窗对应矩形窗。[5]

*理解窗函数和频谱泄露（？）
“窗的波形图显示了窗本身为一个连续的频谱，有一个主瓣，若干旁瓣。 主瓣是时域信号频率成分的中央，旁瓣接近于0。 旁瓣的高度显示了加窗函数对于主瓣周围频率的影响。 对强正弦信号的旁瓣响应可能会超过对较近的弱正弦信号主瓣响应（即，频谱泄漏）。 一般而言，低旁瓣会减少 FFT 的泄漏，但是增加主瓣的带宽。 旁瓣的跌落速率是旁瓣峰值的渐进衰减速率。 增加旁瓣的跌落速率，可减少频谱泄漏。”

“Hamming 窗和 Hanning 窗都有正弦波的外形。 两个窗都会产生宽波峰低旁瓣的结果。 Hanning 窗在窗口的两端都为0，杜绝了所有不连续性。 Hamming 窗的窗口两端不为0，信号中仍然会呈现不连续性。 Hamming 窗擅长减少最近的旁瓣，但是不擅长减少其他旁瓣。 Hamming 窗和 Hanning 适用于对频率精度要求较高对旁瓣要求较低的噪声测量。 ”

二、Matlab 生成信号代码

2.1 信号类型

Create a signal consisting of three noisy sinusoids and a chirp, sampled at 200 kHz for 0.1 second. The frequencies of the sinusoids are 1 kHz, 10 kHz, and 20 kHz. The sinusoids have different amplitudes and noise levels. The noiseless chirp has a frequency that starts at 20 kHz and increases linearly to 30 kHz during the sampling. [1]
Note: 利用矩阵乘法给不同频率成分的信号加权重，代码简单清晰。
Fs = 200e3;
Fc = [1 10 20]'1e3;
Ns = 0.1Fs;

t = (0:Ns-1)/Fs;
x = [1 1/10 10]sin(2piFct)+[1/200 1/2000 1/20]*randn(3,Ns);
x = x+chirp(t,20e3,t(end),30e3);

2.2 数据类型（for SignalAnalyzer APP）

2.2.1 表

duration 数组

timetable 数组

table2timetable
array2timetable = array2table + table2timetable

timeseries 对象

labeledSignalSet 对象

stackedplot（具有公共 x 轴的几个变量的堆叠图）R2023a

2.2.2 组

三、SignalAnalyzer APP

3.1 SignalAnalyzer APP workflow: [4]

选择要分析的信号 - 选择 MATLAB® 工作区中可用的任何信号。该 App 接受具有固有时间信息的数值数组和信号，例如 MATLAB timetable 数组、timeseries 对象和 labeledSignalSet 对象。有关详细信息，请参阅Data Types Supported by Signal Analyzer。使用采样率、数值向量、duration 数组或 MATLAB 表达式向信号添加时间信息。

Preprocess Signals - 使用截断、削波或裁剪操作来编辑信号。低通、高通、带通或带阻滤波器信号。去趋势并计算信号包络。使用移动平均值、回归、Savitzky-Golay 滤波器或其他方法对信号进行平滑处理。使用小波对信号进行去噪。更改信号的采样率或将非均匀采样的信号插值到均匀网格上（重采样）。使用您自己的自定义函数预处理信号。生成 MATLAB 函数来自动执行预处理操作。

探查信号 - 绘制、测量和比较数据、其频谱、频谱图或尺度图。寻找时域、频域和时频域中的特性和模式。计算持久频谱以分析偶发信号，并使用重排来锐化频谱图估计。从信号中提取感兴趣的区域。测量时域中的信号统计量（如最小值、最大值、均值和均方根值）。

共享分析 - 将显示内容作为图像从 App 复制到剪贴板。将信号导出到 MATLAB 工作区或将其保存到 MAT 文件。生成 MATLAB 脚本，以自动计算功率谱、频谱图或持久频谱估计，并自动提取感兴趣的区域。保存信号分析器会话，以便以后或在另一台机器上继续分析。

一个整理得很美观得教程网页 [7]

3.2 COMPARE spectrogram.m & pspectrum(‘spectrogram’)

（两种函数画出来的频谱图的不同之处）
选择好合适的 beta | leakage, 使 Kaiser window ~ hann/hanning window，给出合适的 FrequencyResolution，pspectrum 可以得到同 spectrogram 类似的结果；反之亦然
在 pspectrum 中指定 TimeResolution 没有得到想要的结果（pspectrum 中显示的时间分辨率和 overlap 变化无关，更接近 overlap 为零的情况。所以如果指定时间分辨率，也应该直接用 window/fs 的结果，而非 (window-noverlap)/fs），为什么？pspectrum 中指定频率和时间分辨率的操作有什么不同？Signal Analyzer 中没有改变频率分辨率的选项怎么办？–可以看一下帮助文档[5]中的 RWB 的描述和更细节的处理，暂时搁置@2022/8/26。
类似的分辨率和窗函数，为什么得到的功率值不同？spectrogram 中显示的是 dB/Hz(‘spectrumtype’ 默认为 ‘psd’-功率谱密度, 也可设置为 ‘power’-功率)，pspectrum 中是 dB。
3 ways to sharpen time-frequency estimates
the reassigned spectrogram: pspectrum(“reassign”)
the synchrosqueezed transform ：fsst ifsst
the cross-spectrogram: xspectrogram
COMPARE fsst & pspectrum(‘Reassign’,true)/spectrogram(‘reassigned’) （同步压缩 和 重新分配）
fsst 可进行逆变换来重建模式 - Unlike the reassigned spectrogram, the synchrosqueezed transform is invertible and thus can reconstruct the individual modes that compose the signal.
fsst 只能调节窗口长度和类型（window），不能调节重叠和DFT点数（noverlap 和 nfft）。此限制是为了能够进行逆变换。[10]
DFT点数等于指定窗口的长度。
相邻窗口段之间的重叠比窗口长度少1。
重新分配仅在频率上执行。
spectrogram(sig,kaiser(256,10),255,256,fs);

3.3 COMPARE persistence spectrum & spectrogram （持久频谱 和 频谱图）

persistence spectrum 只存在于 SignalAnalyzer 或等价函数 pspectrum 中，persistence spectrum 与 spectrogram 都是 STFT 的图像显示, 只是图像显示方式不同.
STFT对每一帧的时间序列进行 fft.
spectrogram 用热力图来表达, 横坐标表示该帧对应的时间, 纵坐标表示频率, 颜色表示该帧,该频率下的power(dB).（缺点: 1 时间分辨率, 频率分辨率受"测不准定理"的限制. 2 很多人对相近颜色的区别不是太敏感.）
persistence spectrum 可以理解为对每一帧的频谱线进行了重叠, 即很多个频谱画在了一张图上, 然后对图形进行离散化, 从而可以看出哪些部位重叠度比较高. 最后也是热力图.横坐标表示频率, 纵坐标表示power (dB), 颜色表示重叠的比例(0-100%). 虽然是热力图, 但是可以清晰的看到一条一条"频谱线", 特别是当帧的数量不是很多时.（优点：可以清晰看到由哪些局部信号构成,比靠颜色来分辨 power(dB)直观的多；不足：却无法得知这些"频谱线"的先后顺序.）
spectrogram(STFT) & scalogram(cwt) （频谱图 和 尺度图）
cwt 结果：高频部分低 Fres 高 Tres，低频部分高 Fres 低 Tres
暂时搁置@2022/8/26

四、Matlab 信号预处理

五、Matlab 相关分析

六、Matlab 谱分析

总表和简单介绍

相关、频谱、互谱、功率谱

FUNCTION	FUNCTION of FUNCTION	COMMENT
fft	快速傅里叶变换 (Fast Fourier transform)	频谱
xcorr	互相关 (Cross-correlation)
mscohere	(Magnitude-squared coherence)	? Cxy
csd	互谱密度 (Cross Spectral Density)	csd has been replaced by CPSD.
psd	功率谱密度 (Powder Spectral Density)	psd has been deprecated, use PERIODOGRAM or PWELCH instead.)
cpsd	互功率谱密度 (Cross power spectral density)	CPSD does not support detrending. Please use the DETREND function if you need to detrend your signal.
pwelch	Welch 功率谱密度估计 (Welch’s power spectral density estimate)

注意和统计概念进行区分

FUNCTION	FUNCTION of FUNCTION	COMMENT
xcov	互协方差 ()
corr		Pearson 相关系数 = 两个变量的协方差除以标准差的乘积。3

时频谱

非平稳信号是频率成分随时间变化的信号。
非平稳信号的频谱图是对其频率成分的时间演变的估计，应该使用快速傅里叶变换(STFT) 或者是连续小波变换得到。
STFT 相关函数：
spectrogram: Spectrogram using short-time Fourier transform
xspectrogram: Cross-spectrogram using short-time Fourier transforms 输出结果只包含互相关谱的强度信息，不包含相位信息；需要相位信息需要自己编写程序
pspectrum: Analyze signals in the frequency and time-frequency domains 有不同“模式”可供选择，是信号分析(SignalAnalyzer) APP 的等价函数
fsst: Fourier synchrosqueezed transform 傅里叶同步压缩变换？
tfridge: Time-frequency ridges

谱估计相关函数

其他函数

pow2db: Convert power to decibels
downsample: Decrease sample rate by integer factor

常见输入变量的设定

计算相关

window
noverlap
nfft

加窗前分辨率的计算公式：
fs 采样率（1s 采集点数）
N 一次输入 FFT 的采样点数，（对应 nfft）
T = N/fs 采样时长（这里默认窗口覆盖选取的全部时间）
dt = 1/fs 时间分辨率（t = (0:N-1)dt = 0:(1/fs):(T-1/fs)）
df = 1/T = 1/(Ndt) = fs/N 频率分辨率（f = (0:N-1)df = (0:N-1)(fs/N) = 0:(1/T):(fs-1/T)），选择频率段 0:fre 对应索引 1: fre/df = 1:freT
dtdf = 1/N
加窗后分辨率的计算公式[9]：
df = fs/nfft
对应 [nfft/2] 或 [nfft/2+1] frequencies
dt = (t(end)-t(1))/((length(x)-noverlap)/(window-noverlap)) = (window-noverlap)/fs
对应 [(length(x)-noverlap)/(window-noverlap)] time bins
注意：上面的 t 是频谱图的输出，不是原始的 t
参数选择总结：
影响频率分辨率：nfft，window 形状/ leakage，补零（？）
影响时间分辨率：window width，noverlap,选择的窗口内近似为平稳信号即可，window-noverlap 影响最后呈现效果
% fs = 1e6
% nfft = 10248
% nT = nfft/fs
% window = 10248
% noverlap = round(window*0.9)
%
% df = fs/nfft
% dt = 1/fs
% dwt = (window - noverlap)/fs
% dwtno = window/fs

显示相关

w
f
fs
freqrange
trace    'mean'（默认）'minhold' or 'maxhold' 分别对应各片段(segments) 功率谱估计的的平均值、最小值和最大值

本人原笔记：谱估计方法及 MATLAB 实现
参考资料
[1]C:\Users\User\Documents\MATLAB\Examples\R2020b\signal\WelchSpectrumEstimatesExample\WelchSpectrumEstimatesExample.mlx
[2] MATLAB数字信号处理（1）四种经典功率谱估计方法比较_FPGADesigner的博客-CSDN博客_数字信号功率谱
[3] Matlab 实现经典功率谱分析和估计_MissXy_的博客-CSDN博客_matlab 功率谱
[4] 使用信号分析器 - MATLAB & Simulink - MathWorks 中国
[5] pspectrum matlab 帮助文档
[6] 信号分析器中的频谱图计算 - MATLAB & Simulink - MathWorks 中国
[7] Matlab Signal Analyzer 信号分析 | 萤火 (jychen.cn) - 似乎要教育网才能打开
[8] 信号频域分析方法的理解（频谱、能量谱、功率谱、倒频谱、小波分析） - 知乎 (zhihu.com)
[9]spectrogram matlab 帮助文档
[10] fsst matlab 帮助文档

---

1. Matlab Signal Analyzer 信号分析 | 萤火 (jychen.cn) - 似乎要教育网才能打开 ↩︎

2. 频域特征值提取的MATLAB代码实现（频谱、功率谱、倒频谱） ↩︎

3. 深入理解Pearson相关系数和Matlab corr()、corrcoef()仿真 ↩︎