Syslab使用手册
⋮
2026b

# xspectrogram

使用短时傅立叶变换的互谱图

函数库: TySignalProcessing

提示

推荐使用 ty_xspectrogram 以获取更好的性能表现和使用体验。

# 语法

s, = xspectrogram(x, y)
s, = xspectrogram(x, y, window)
s, = xspectrogram(x, y, window, noverlap)
s, = xspectrogram(x, y, window, noverlap, nfft)
s, w, t = xspectrogram(___)
s, f, t = xspectrogram(___, fs)
s, w, t = xspectrogram(x, y, window, noverlap, w)
s, f, t = xspectrogram(x, y, window, noverlap, f, fs)
___, c = xspectrogram(___)
___ = xspectrogram(___, freqrange)
___ = xspectrogram(___, Name, Value)
___ = xspectrogram(___, spectrumtype)
xspectrogram(___)
xspectrogram(___, freqloc)

# 说明

s, = xspectrogram(x, y) 返回 x 和 y 指定的信号的互谱图。输入信号必须是元素数量相同的向量。s 的每一列都包含 x 和 y 共同的短期、时间局部化频率内容的估计。

s, = xspectrogram(x, y, window) 使用 window 将 x 和 y 划分为多个分段并执行加窗操作。

s, = xspectrogram(x, y, window, noverlap) 使用相邻分段之间重叠的 noverlap 样本。

s, = xspectrogram(x, y, window, noverlap, nfft) 使用 nfft 采样点来计算离散傅立叶变换。

s, w), t = xspectrogram(___) 返回归一化频率的向量w和计算互谱图的时刻的向量 t。此语法可以包括以前语法中输入参数的任何组合。

s, f, t = xspectrogram(___, fs) 返回频率向量 f，用采样率 fs 表示。fs 必须是 xsspectrogram 的第六个输入。要输入采样率并仍然使用前面可选参数的默认值，请将这些参数指定为空。

s, w), t = xspectrogram(x, y, window, noverlap, w) 返回 w 中指定的归一化频率下的互谱图。

s, f, t = xspectrogram(x, y, window, noverlap, f, fs) 返回 f 中指定频率的互谱图。

___, c = xspectrogram(___) 还返回一个矩阵 c，该矩阵包含输入信号的时变复互谱的估计。互谱图 s 是 c 的幅值。

___ = xspectrogram(___, freqrange) 返回 freqrange 指定频率范围内的互谱图。freqrange 的有效选项为 "onesided"、"twosided" 和 "centered"。

___ = xspectrogram(___, Name, Value) 使用名称-值对参数指定其他选项。选项包括最小阈值和输出时间维度。

___ = xspectrogram(___, spectrumtype) 如果频谱类型指定为 "psd"，则返回短期互功率谱密度估计值；如果频谱类型指定为 "power"，则返回短期互功率谱估计值。

xspectrogram(___) 在没有输出参数的情况下，在当前图形窗中绘制互谱图。

xspectrogram(___, freqloc) 指定绘制频率的轴。将 freqloc 指定为 "xaxis" 或 "yaxis"。

# 示例

线性线性调频互谱图

生成以 1 MHz 采样 10 毫秒的两个线性线性调频。

  • 第一线性调频具有 150 kHz 的初始频率，该初始频率在测量结束时增加到 350 kHz；

  • 第二线性调频具有 200 kHz 的初始频率，该初始频率在测量结束时增加到 300 kHz。

添加高斯白噪声，使信噪比为 40 dB。

using TyMath
using TySignalProcessing
using TyPlot
using TyControlSystems
nSamp = 10000
Fs = 1000e3
SNR = 40
t = [0:(nSamp - 1);] ./ Fs

rng = MT19937ar(1234)
x1 = chirp(t, 150e3, t[end], 350e3)
x1 = x1 + randn(rng, size(x1)) * std(x1) / db2mag(SNR)

x2 = chirp(t, 200e3, t[end], 300e3)
x2 = x2 + randn(rng, size(x2)) * std(x2) / db2mag(SNR)

计算并绘制两个线性调频信号的互谱图。将信号划分为 200 个采样段，并用汉明窗对每个信号段进行加窗处理。指定相邻段之间有 80 个样本重叠，DFT 长度为 1024 个样本。

figure()
xspectrogram(x1, x2, hamming(200), 80, 1024, Fs, "yaxis"; plotfig=true)

修改第二个线性调频信号，使频率在测量过程中从 50 kHz 上升到 350 kHz。使用形状因子 β=5 的 500 样本 Kaiser 窗来对分段加窗。指定 450 个重叠样本和 256 的 DFT 长度。计算并绘制互谱图。

x2 = chirp(t, 50e3, t[end], 350e3)
x2 = x2 + randn(rng, size(x2)) * std(x2) / db2mag(SNR)

figure()
xspectrogram(x1, x2, kaiser(500, 5), 450, 256, Fs, "yaxis"; plotfig=true)

在这两种情况下，该函数都能突出显示两个信号的共同频率内容。

语音信号互谱图

加载一个包含两个语音信号的文件，采样频率为 44 100 Hz。

  • 第一个信号是一个女声说 "transform function" 的录音；

  • 第二个信号是同一个女声说 "reform justice." 的录音。

绘制这两个信号的曲线图。将第二个信号垂直偏移，使两个信号都可见。

using TySignalProcessing
using TyBase
using TyPlot
pkg_dir = pkgdir(TySignalProcessing)
source_path = pkg_dir * "/examples/Resource/voice.mat"
load(source_path)
# To hear, type using TyDSPSystem;soundsc(transform,fs);pause(2);soundsc(reform,fs)
t = [0:(length(reform) - 1);] / fs

plot(t, transform, t, reform .+ 0.3)
legend(["\"Transform function\"", "\"Reform justice\""])

计算两个信号的互谱图。将信号分成 1000 个采样段，并用汉明窗对其进行加窗处理。在相邻信号段之间指定 800 个重叠样本。仅包含 4 kHz 以下的频率。

nwin = 1000
nvlp = 800
fint = 0:4000

figure()
s, f, t = xspectrogram(
    transform, reform, hamming(nwin), nvlp, fint, fs, "yaxis"; plotfig=true
)

互谱图突出显示了信号共同频率较高的时间间隔。音节 "form" 尤其明显。

两个二次线性调频之间的相移

产生两个二次线性调频，每个以 1 kHz 采样，持续 2 秒钟。两个线性调频的初始频率都是 100 赫兹，在测量中途增加到 200 赫兹。与第一个线性调频声相比，第二个线性调频声的相位差为 23°。

using TySignalProcessing
using TyPlot
fs = 1e3
t = 0:(1 / fs):2

y1 = chirp(t, 100, 1, 200, "quadratic", 0)
y2 = chirp(t, 100, 1, 200, "quadratic", 23)

计算线性调频信号的复交谱图，提取它们之间的相移。将信号分成 128 个采样段。在相邻段之间指定 120 个重叠样本。使用形状因子 β = 18 的 Kaiser 窗为每个分段加窗，并指定 128 个样本的 DFT 长度。使用 xspectrogram 的绘图功能显示互谱图。

_, f, xt, c = xspectrogram(y1, y2, kaiser(128, 18), 120, 128, fs)

figure()
xspectrogram(y1, y2, kaiser(128, 18), 120, 128, fs, "yaxis"; plotfig=true)
复信号互谱图

产生两个信号，每个信号的采样频率为 3 kHz，持续 1 秒钟。第一个信号是二次线性调频信号，其频率在测量过程中从 300 Hz 增至 1300 Hz。该线性调频信号包含在白高斯噪声中。第二个信号也包含在白噪声中，是一个频率内容正弦变化的线性调频信号。

using TyMath
using TySignalProcessing
using TyPlot
fs = 3000
t = 0:(1 / fs):(1 - 1 / fs)
rng = MT19937ar(1234)
x1 = chirp(t, 300, t[end], 1300, "quadratic") + randn(rng, size(t)) / 100

x2 = exp.(2im * pi * 100 * cos.(2 * pi * 2 * t)) + randn(rng, size(t)) / 100

计算并绘制两个信号的互谱图。将信号分成 256 个样本段，相邻信号段之间重叠 255 个样本。使用形状系数 β = 30 的 Kaiser 窗对信号段进行分段。使用默认的 DFT 点数。将横谱图以零频率为中心。

nwin = 256

figure()
xspectrogram(x1, x2, kaiser(nwin, 30), nwin - 1, [], fs, "centered", "yaxis"; plotfig=true)

计算功率谱，而不是功率谱密度。将小于 -40 dB 的值设为零。以奈奎斯特频率为中心绘制。

figure()
xspectrogram(
    x1,
    x2,
    kaiser(nwin, 30),
    nwin - 1,
    [],
    fs,
    "power",
    "MinThreshold",
    -40,
    "yaxis";
    plotfig=true,
)
title("Cross-Spectrogram of Quadratic Chirp and Complex Chirp")

阈值处理进一步突出了共同频率区域。

两个序列的互谱图

计算并绘制两个序列的互谱图。

指定每个序列的长度为 4096 个样本。

using TyMath
using TySignalProcessing
using TyPlot
N = 4096

要创建第一个序列，需要生成一个嵌入白高斯噪声的凸二次线性调频，并对其进行带通滤波。

  • 线性调频的初始归一化频率为 0.1π，在测量结束时增加到 0.8π；

  • 16 阶滤波器通过 0.2π 和 0.4π rad/sample之间的归一化频率，其阻带衰减为 60 dB。

rng = MT19937ar(1234)
rx = chirp(0:(N - 1), 0.1 / 2, N, 0.8 / 2, "quadratic", [], "convex") + randn(rng, N) / 100
b1, a1 = cheby2(16, 60, [0.2 0.4], "bandpass", "z")
x, = filter1(real(b1), a1, rx)

要创建第二个序列，需要生成一个嵌入白高斯噪声的线性线性调频，并对其进行带阻滤波。

  • 线性调频的初始归一化频率为 0.9π，在测量结束时降至 0.1π；

  • 16 阶滤波器可阻止 0.6π 和 0.8π rad/sample之间的归一化频率，通带纹波为 1 dB。

ry = chirp(0:(N - 1), 0.9 / 2, N, 0.1 / 2) + randn(rng, N) / 100
b2, a2 = cheby1(16, 1, [0.6 0.8], "bandstop", "z")
y, = filter1(b2, a2, ry)

绘制两个序列。垂直偏移第二个序列，使两个序列都可见。

plot(x)
hold("on")
plot(y .+ 2)

计算并绘制 x 和 y 的互谱图。指定相邻片段之间有 500 个重叠样本和 2048 个 DFT 点。

figure()
xspectrogram(x, y, hamming(512), 500, 2048, "yaxis"; plotfig=true)

将小于 -50 dB 的互谱图值设为零。

figure()
xspectrogram(x, y, hamming(512), 500, 2048, "MinThreshold", -50, "yaxis"; plotfig=true)

频谱图显示了滤波器增强或抑制的频率区域。

# 输入参数

x - 输入信号
向量

以向量形式指定的输入信号。

示例： cos.(pi/4*(0:159)) .+ randn(160) 表示嵌入白高斯噪声的正弦波。

y - 输入信号
向量

以向量形式指定的输入信号。

示例： cos.(pi/4*(0:159)) .+ randn(160) 表示嵌入白高斯噪声的正弦波。

window - 窗
整数 | 向量

窗，指定为整数或行或列向量。使用窗将信号划分为若干段：

  • 如果 window 为整数，则 spectrogram 会将 x 分割成长度为 window 的段落，并用该长度的汉明窗对每个分段进行加窗处理；

  • 如果 window 是向量，则 spectrogram 将 x 分割成与向量长度相同的段，并用 window 对每个分段加窗。

如果 x 的长度无法精确划分为整数个无重叠采样的分段，则会对 x 进行相应的截断。

如果将 window 指定为空，则 spectrogram 会使用一个 Hamming 窗，将 x 分成 8 个无重叠采样点的分段。

noverlap - 重叠的样本数
正整数

重叠样本数，以正整数表示。

  • 如果 window 是标量，则 noverlap 必须小于 window；

  • 如果 window 是向量，则 noverlap 必须小于 window 的长度。

如果指定 noverlap 为空，则 spectrogram 会使用一个能产生 50% 片段重叠的数字。如果未指定段长度，函数会将 noverlap 设为 ，其中 是输入信号的长度， 符号表示下限函数。

nfft - DFT 点数
正整数标量

DFT 点数，指定为正整数标量。如果指定 nfft 为空，则 spectrogram 会将参数设置为 ，其中 。

  • 如果 window 是标量，则 = window；

  • 如果 window 是一个向量，则 = length(window)。

w - 归一化频率
向量

指定为向量的归一化频率。w 必须至少有两个元素，否则函数会将其解释为 nfft。归一化频率单位为 rad/sample。

示例： pi./[2 4］

f - 频率
向量

频率，指定为向量。f 必须至少有两个元素。f 的单位由采样率 fs 指定。

fs - 采样率
1 Hz（默认） | 正标量

采样率，指定为正标量。采样率是单位时间内的采样次数。如果时间单位是秒，那么采样率的单位就是赫兹。

freqrange - 互谱估计的频率范围
"onesided" | "twosided" | "centered"

互谱估计的频率范围，指定为 "onesided"、"twosided" 或 "centered"。对于实值信号，默认为 "onesided"。对于复值信号，默认为 "twosided"，指定 "onesided" 会导致错误。

  • "onesided" - 返回实数输入信号的单边互谱图。如果 nfft 为偶数，则 s 有 nfft/2 + 1 行，并在 [0, π] rad/sample区间内计算。如果 nfft 为奇数，则 s 有 (nfft + 1)/2 行，区间为 [0, π] rad/sample。如果指定 fs，则间隔分别为 [0, fs/2] 周期/单位时间和 [0, fs/2) 周期/单位时间；

  • "twosided" - 返回实数或复数信号的双侧互谱图。s 有 nfft 行，并且是在区间 [0，2π) rad/sample 上计算的。如果指定 fs，则区间为 [0，fs) 周期/单位时间；

  • "centered" - 返回实数或复数信号的居中双面互谱图。如果 nfft 为偶数，则 s 在 (-π, π] rad/sample 的区间内计算。如果 nfft 为奇数，则 s 在 (-π, π) rad/sample 上计算。如果指定 fs，则间隔时间分别为 (-fs/2, fs/2] 周期/单位时间和 (-fs/2, fs/2) 周期/单位时间。

spectrumtype - 互功率谱缩放
"psd"（默认） | "power"

互功率谱缩放，指定为 "psd "或 "power"。

  • 省略频谱类型或指定 "psd"，将返回互功率谱密度；

  • 指定 "power" 时，互功率谱密度的每个估计值都会按分辨率带宽缩放，分辨率带宽取决于窗的等效噪声带宽和片段持续时间。结果是每个频率的功率估计值。

freqloc - 频率显示轴
"xaxis"（默认） | "yaxis"

频率显示轴，指定为 "xaxis" 或 "yaxis"。

"xaxis" - 在 x 轴上显示频率，在 y 轴上显示时间。

"yaxis" -在 y 轴上显示频率，在 x 轴上显示时间。

如果调用带有输出参数的 xspectrogram，此参数将被忽略。

# 名称-值对参数

指定可选的以逗号分隔的 Name,Value 参数对。Name 是参数名称，Value 是相应的值。您可以按任意顺序指定多个名称和值对参数，如 Name1,Value1,...,NameN,ValueN。

示例： xspectrogram(x,100,"OutputTimeDimension","downrows") 将 x 和 y 分割成长度为 100 的分段，并用该长度的 Hamming 窗对每个分段进行加窗处理。频谱图的输出时间维度为下行。

"MinThreshold" - 阈值
-Inf（默认） | 实标量

阈值，指定为 MinThreshold 和以分贝为单位的实数标量。xspectrogram 将 10 log10(s) ≤ thresh 的 s 元素置零。

"OutputTimeDimension" - 输出时间维度
acrosscolumns（默认） | downrows

输出时间维度，指定为由 OutputTimeDimension 和 acrosscolumns 或 downrows 组成的逗号分隔对。如果希望 s、ps、fc 和 tc 的时间维度沿行显示，频率维度沿列显示，则将此值设为 downrows。如果希望 s、ps、fc 和 tc 的时间维度跨列，频率维度沿行，则将此值设为 acrosscolumns。如果调用函数时没有输出参数，则忽略此输入。

# 输出参数

s - 互谱图
矩阵

互谱图，以矩阵形式返回。时间在 s 的各列之间递增，频率从 0 开始在各行之间递增。

  • 如果输入信号 x 和 y 的长度为 N，则 s 有 k 列，其中：

    如果 window 是标量，则 k = ⌊（N - noverlap）/（window - noverlap）⌋；

    如果 window 是向量，则 k = ⌊（N - noverlap）/（length(window) - noverlap）⌋。

  • 如果输入信号为实数，且 nfft 为偶数，则 s 有 (nfft/2 + 1) 行；

  • 如果输入信号为实数，且 nfft 为奇数，则 s 有 (nfft + 1)/2 行；

  • 如果输入信号是复数，则 s 有 nfft 行。

w - 归一化频率
向量

归一化频率，以向量形式返回。w 的长度等于 s 的行数。

t - 时间瞬时
向量

时间瞬时，以向量形式返回。t 中的时间值对应于使用 window 指定的每个片段的中点。

f - 周期频率
向量

周期频率，以向量形式返回。f 的长度等于 s 的行数。

c - 时变复数互功率谱
矩阵

时变复数互功率谱，以矩阵形式返回。互谱图 s 是 c 的幅度。

# 参考文献

[1] Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.

[2] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

# 另请参阅

cpsd | mscohere | spectrogram