# ty_spectrogram
使用短时傅立叶变换的频谱图
函数库: TySignalProcessing
# 语法
s, = ty_spectrogram(x)
s, = ty_spectrogram(x, window)
s, = ty_spectrogram(x, window, noverlap)
s, = ty_spectrogram(x, window, noverlap, nfft)
s, w, t = ty_spectrogram(___)
s, f, t = ty_spectrogram(___, fs)
s, w, t = ty_spectrogram(x, window, noverlap, w)
s, f, t = ty_spectrogram(x, window, noverlap, f, fs)
_, _, _, ps = ty_spectrogram(___)
_ = ty_spectrogram(___; reassign = true)
_, _, _, ps, fc, tc = ty_spectrogram(___; reassign = true)
___ = ty_spectrogram(___; freqrange = String)
___ = ty_spectrogram(___; Name = Value)
___ = ty_spectrogram(___stype = String)
ty_spectrogram(___; plotfig = true)
ty_spectrogram(___; plotfig = true, freqloc = String)
# 说明
s, = ty_spectrogram(x) 返回输入信号 x 的短时傅立叶变换（STFT）。s 的每列包含 x 的短期、时间局部化频率内容的估计。s 的幅度平方被称为 x[1] 的频谱图时间-频率表示。
s, = ty_spectrogram(x, window) 使用 window 将信号划分为数段并加窗。
s, = ty_spectrogram(x, window, noverlap) 使用相邻段之间重叠的 noverlap 样本。
s, = ty_spectrogram(x, window, noverlap, nfft) 使用 nfft 采样点计算离散傅立叶变换。
s, w, t = ty_spectrogram(___) 返回归一化频率向量 w 和计算 STFT 的时间时刻向量 t。该语法可以包含前面语法中输入参数的任意组合。
s, f, t = ty_spectrogram(___, fs) 返回频率的向量 f，用采样率 fs 表示。fs 必须是 ty_spectorgram 的第五个输入。要输入采样率并仍然使用前面可选参数的默认值，请将这些参数指定为空，[]。
s, w, t = ty_spectrogram(x, window, noverlap, w) 以 w 中指定的归一化频率返回 STFT。w 必须至少有两个元素，否则函数会将其解析为 nfft。
s, f, t = ty_spectrogram(x, window, noverlap, f, fs) 以 f 中指定的频率返回 STFT。f 必须至少有两个元素，否则函数会将其解析为 nfft。
_, _, _, ps = ty_spectrogram(__) 还会返回一个与 x 的频谱图成正比的矩阵 ps。
如果指定频谱类型为 "psd"，则 ps 的每一列都包含一个窗分段的功率谱密度（PSD）估计值；
如果指定频谱类型为 "power"，则 ps 的每一列都包含一个窗分段的功率谱估计值。
_ = ty_spectrogram(___; reassign = true) 将每个 PSD 或功率谱估算值重新分配到其能量中心的位置。如果您的信号包含定位良好的时间或频谱成分，那么该选项生成的频谱图会更清晰。
_, _, _, ps, fc, tc = ty_spectrogram(__) 还会返回两个矩阵 fc 和 tc，其中包含每个 PSD 或功率谱估计值的能量中心频率和时间。
___ = ty_spectrogram(___; freqrange = String) 返回 freqrange 指定频率范围内的 PSD 或功率谱估计值。freqrange 的有效选项为 "onesided"、"twosided" 和 "centered"。
___ = ty_spectrogram(___; Name = Value) 使用名称值参数指定其他选项。选项包括最小阈值和输出时间维度。
___ = ty_spectrogram(___; stype = String) 如果频谱类型指定为 "psd"，则返回短期互功率谱密度估计值；如果频谱类型指定为 "power"，则返回短期互功率谱估计值。
ty_spectrogram(___; plotfig = true) 在输入 plotfig = true 的情况下，在当前图形窗中绘制互功率谱图。
ty_spectrogram(___; plotfig = true, freqloc = String) 指定绘制频率的轴。将 freqloc 指定为 "xaxis" 或 "yaxis"。
# 示例
频谱图的默认值
生成
using TySignalProcessing
using TyPlot
using TyMath
N = 1024
n = 0:(N - 1)
w0 = 2 * pi / 5
x = @. sin(w0 * n) + 10 * sin(2 * w0 * n)
使用默认函数计算短时傅立叶变换。绘制频谱图。
s, = ty_spectrogram(x)
figure()
ty_spectrogram(x; plotfig=true)
重复计算。
将信号分成长度为
的段落； 使用 Hamming 窗对各部分进行窗处理；
指定连续部分之间的重叠率为 50%；
计算 FFT 时，使用
点，其中 。
验证这两种方法得到的结果是否相同。
Nx = length(x)
nsc = floor(Int, Nx / 4.5)
nov = floor(Int, nsc / 2)
nff = max(256, 2^nextpow2(nsc))
t, = ty_spectrogram(x, hamming(nsc), nov, nff)
maxerr = maximum(@. abs(abs(t) - abs(s)))
maxerr = 0.0
将信号分成 8 个长度相等的部分，各部分之间重叠 50%。指定与上一步相同的 FFT 长度。计算短时傅里叶变换，并验证其结果是否与前两个步骤相同。
ns = 8
ov = 0.5
lsc = floor(Int, Nx / (ns - (ns - 1) * ov))
t, = ty_spectrogram(x, lsc, floor(Int, ov * lsc), nff)
maxerr = maximum(@. abs(abs(t) - abs(s)))
maxerr = 0.0
沿 x 轴的频率
产生一个二次线性调频信号 x，以 1 kHz 的频率采样，持续 2 秒钟。线性调频信号的频率最初为 100 赫兹，在 t = 1 秒时越过 200 赫兹。
using TySignalProcessing
using TyPlot
t = 0:0.001:2
x = chirp(t, 100, 1, 200, "quadratic")
计算并显示 x 的频谱图。
将信号划分为长度为 128 的段落，并使用汉明窗；
在相邻部分之间指定 120 个重叠样本；
在
个频率和 个时间分段评估频谱。
figure()
ty_spectrogram(x, 128, 120, 128, 1e3; plotfig=true, freqloc="xaxis")
用 Blackman 窗取代汉明窗。将重叠减少到 60 个样本。绘制时间轴，使其数值从上到下递增。
figure()
ty_spectrogram(x, blackman(128), 60, 128, 1e3; plotfig=true, freqloc="xaxis")
Ylim = ylim()
ylim([maximum(Ylim), minimum(Ylim)])
线性调频的功率谱密度
计算并显示从 100 Hz 开始、在 t = 1 秒时越过 200 Hz 的二次线性调频的每个分段的 PSD。指定采样率为 1 kHz，段长度为 128 个采样点，重叠长度为 120 个采样点。使用 128 个 DFT 点和默认的 Hamming 窗。
using TySignalProcessing
using TyPlot
fs = 1000
t = 0:(1 / fs):2
x = chirp(t, 100, 1, 200, "quadratic")
figure()
ty_spectrogram(x, 128, 120, 128, fs; plotfig=true)
title("Quadratic Chirp")
计算并显示以 1 kHz 频率采样的线性线性调频信号（从直流开始，在 t = 1 秒时越过 150 Hz）各分段的 PSD。指定段长度为 256 个采样点，重叠长度为 250 个采样点。使用默认的 Hamming 窗和 256 个 DFT 点。
x = chirp(t, 0, 1, 150)
figure()
ty_spectrogram(x, 256, 250, 256, fs; plotfig=true)
title("Linear Chirp")
计算并显示以 1 kHz 频率采样的对数线性调频信号（从 20 Hz 开始，在 t = 1 秒时越过 60 Hz）各分段的 PSD。指定段长度为 256 个采样点，重叠长度为 250 个采样点。使用默认的 Hamming 窗和 256 个 DFT 点。
x = chirp(t, 20, 1, 60, "logarithmic")
figure()
ty_spectrogram(x, 256, 250, [], fs; plotfig=true)
title("Logarithmic Chirp")
频率轴使用对数刻度。频谱图变成了一条直线。
gca().set_yscale("log")
gca().set_ylim([4, 500])
频谱图和瞬时频率
使用频谱图功能测量和跟踪信号的瞬时频率。
生成一个采样频率为 1 kHz 的二次线性调频信号，持续两秒钟。指定线性调频信号的频率，使其最初为 100 赫兹，一秒后增加到 200 赫兹。
using TySignalProcessing
using TyPlot
fs = 1000
t = 0:(1 / fs):(2 - 1 / fs)
y = chirp(t, 100, 1, 200, "quadratic")
使用频谱图功能中的短时傅里叶变换估算线性调频信号的频谱。将信号划分为长度为 100 的部分，并使用汉明窗。在相邻部分之间指定 80 个重叠样本，并评估 [100/2 + 1]=51 频率的频谱。
figure()
ty_spectrogram(y, 100, 80, 100, fs; plotfig=true)
复信号频谱图
生成 512 个正弦变化频率的线性调频信号样本。
using TySignalProcessing
using TyPlot
N = 512
n = 0:(N - 1)
x = @. exp(1im * pi * sin(8 * n / N) * 32)
计算线性调频信号的居中两边短时傅里叶变换。将信号分为 32 个采样段，其中 16 个采样段重叠。指定 64 个 DFT 点。绘制频谱图。
scalar, fs, ts = ty_spectrogram(x, 32, 16, 64; freqrange="centered")
figure()
ty_spectrogram(x, 32, 16, 64; freqrange="centered", plotfig=true)
在 (-π,π] 间隔内计算 64 个等距频率的频谱图，也能得到相同的结果。无需使用 "centered" 选项。
fintv = (-pi + pi / 32):(pi / 32):pi
vector, fv, tv = ty_spectrogram(x, 32, 16, fintv)
figure()
ty_spectrogram(x, 32, 16, fintv; plotfig=true)
二次线性调频的重新分配频谱图
生成一个线性调频信号，以 1 kHz 的频率采样 2 秒钟。指定线性调频信号的频率，使其最初为 100 赫兹，1 秒后增加到 200 赫兹。
using TySignalProcessing
using TyPlot
Fs = 1000
t = 0:(1 / Fs):2
y = chirp(t, 100, 1, 200, "quadratic")
估计信号的重配频谱图。
将信号划分为长度为 128 的段落，使用形状参数为 β=18 的 Kaiser 窗进行窗化；
在相邻部分之间指定 120 个重叠样本；
在 [128 / 2]=65 个频率和 [(length(x)-120) / (128-120)]=235 个时间分段评估频谱。
figure()
ty_spectrogram(y, kaiser(128, 18), 120, 128, Fs; reassign=true, plotfig=true)
频谱图重新分配和阈值处理
生成一个 1024 Hz 采样信号，持续 2 秒钟。
using TySignalProcessing
using TyPlot
using TyMath
using TyControlSystems
nSamp = 2048
Fs = 1024
t = (0:(nSamp - 1)) / Fs
在第一秒，信号由一个 400 Hz 正弦波和一个凹二次线性调频组成。指定线性调频信号，使其围绕间隔中点对称，起始和终止频率均为 250 赫兹，最低频率为 150 赫兹。
t1 = t[1:Int(nSamp / 2)]
x11 = @. sin(2 * pi * 400 * t1)
x12 = chirp(t1 .- t1[Int(nSamp / 4)], 150, nSamp / Fs, 1750, "quadratic")
x1 = x11 + x12
信号的其余部分由两个频率递减的线性线性调频组成。其中一个线性调频的初始频率为 250 赫兹，随后下降到 100 赫兹。另一个线性调频的初始频率为 400 赫兹，随后降至 250 赫兹。
t2 = t[Int(nSamp / 2 + 1):nSamp]
x21 = chirp(t2, 400, nSamp / Fs, 100)
x22 = chirp(t2, 550, nSamp / Fs, 250)
x2 = x21 + x22
为信号添加白高斯噪声。指定信噪比为 20 dB。重置随机数生成器，以获得可重复的结果。
SNR = 20
rng = MT19937ar(1234)
sig = [x1; x2]
sig = sig + randn(rng, size(sig)) * std(sig) / db2mag(SNR)
计算并绘制信号频谱图。指定一个长度为 63、形状参数为 β=17 的 Kaiser 窗，在相邻部分之间减少 10 个重叠样本，FFT 长度为 256。
nwin = 63
wind = kaiser(nwin, 17)
nlap = nwin - 10
nfft = 256
figure()
ty_spectrogram(sig, wind, nlap, nfft, Fs; plotfig=true)
对频谱图进行阈值处理，将数值小于信噪比的元素置零。
figure()
ty_spectrogram(sig, wind, nlap, nfft, Fs; MinThreshold=-SNR, plotfig=true)
将每个 PSD 估计值重新分配到其能量中心位置。
figure()
ty_spectrogram(sig, wind, nlap, nfft, Fs; reassign=true, plotfig=true)
对重新分配的频谱图进行阈值处理，将值小于信噪比的元素设置为零。
figure()
ty_spectrogram(sig, wind, nlap, nfft, Fs; reassign=true, MinThreshold=-SNR, plotfig=true)
跟踪音频信号中的线性调频
加载一个音频信号，其中包含两个递减的线性调频声和一个宽带飞溅声。计算短时傅立叶变换。将波形划分为 400 个采样段，重叠 300 个采样段。绘制频谱图。
using TySignalProcessing
using TyPlot
using TyBase
pkg_dir = pkgdir(TySignalProcessing)
source_path = pkg_dir * "/examples/Resource/splat.mat"
load(source_path)
# To hear, type using TyDSPSystem;soundsc(y,Fs)
sg = 400
ov = 300
figure()
ty_spectrogram(y, sg, ov, [], Fs; plotfig=true)
使用频谱图功能输出信号的功率谱密度 (PSD)。
s, f, t, p = ty_spectrogram(y, sg, ov, [], Fs)
使用 medfreq 功能跟踪这两个鸣声。要找到较强的低频鸣叫，搜索范围应限制在 100 Hz 以上的频率和宽带声音开始前的时间。
f1 = f .> 100
t1 = t .< 0.75
m1, = medfreq(p[f1, vec(t1)], f[f1])
要找到微弱的高频鸣叫，搜索范围应限制在 2500 赫兹以上，时间在 0.3 秒到 0.65 秒之间。
f2 = f .> 2500
t2 = t .> 0.3 .&& t .< 0.65
m2, = medfreq(p[f2, vec(t2)], f[f2])
将结果叠加到频谱图上。将频率值除以 1000，以千赫为单位表示。
hold("on")
plot(t[t1], m1 / 1000; linewidth=4)
plot(t[t2], m2 / 1000; linewidth=4)
hold("off")
三维频谱图可视化
生成两秒钟的 10 kHz 采样信号。指定信号的瞬时频率为时间的三角函数。
using TySignalProcessing
using TyPlot
fs = 10e3
t = 0:(1 / fs):2
x1 = vco(sawtooth(2 * pi * t, 0.5), [0.1 0.4] * fs, fs)
计算并绘制信号的频谱图。使用长度为 256、形状参数 β=5 的 Kaiser 窗。指定 220 个节间重叠样本和 512 个 DFT 点。在 Y 轴上绘制频率图。使用默认色图和视图。
figure()
ty_spectrogram(x1, kaiser(256, 5), 220, 512, fs; plotfig=true);
# 输入参数
x - 输入信号向量
输入信号，指定为行或列向量。
示例： cos.(pi/4*(0:159)) .+ randn(1,160) 表示嵌入白高斯噪声的正弦波。
window - 窗整数 | 向量
窗，指定为整数或行或列向量。使用窗将信号划分为若干段：
如果 window 为整数，则 ty_spectrogram 会将 x 分割成长度为 window 的段落，并用该长度的汉明窗对每个分段进行窗处理；
如果 window 是向量，则 ty_spectrogram 将 x 分割成与向量长度相同的段，并用 window 对每个分段加窗。
如果 x 的长度无法精确划分为整数个无重叠采样的分段，则会对 x 进行相应的截断。
如果将 window 指定为空，则 ty_spectrogram 会使用一个 Hamming 窗，将 x 分成 8 个无重叠采样点的分段。
noverlap - 重叠的样本数正整数
重叠样本数，以正整数表示。
如果 window 是标量，则 noverlap 必须小于 window；
如果 window 是向量，则 noverlap 必须小于 window 的长度。
如果指定 noverlap 为空，则 ty_spectrogram 会使用一个能产生 50% 片段重叠的数字。如果未指定段长度，函数会将 noverlap 设为
nfft - DFT 点数正整数标量
DFT 点数，指定为正整数标量。如果指定 nfft 为空，则 ty_spectrogram 会将参数设置为
如果 window 是标量，则
= window； 如果 window 是一个向量，则
= length(window)。
w - 归一化频率向量
指定为向量的归一化频率。w 必须至少有两个元素，否则函数会将其解析为 nfft。归一化频率单位为 rad/sample。
示例： pi./[2 4］
f - 周期频率向量
f 必须至少有两个元素，否则函数会将其解析为 nfft。f 的单位由采样率 fs 指定。
fs - 采样率1 Hz（默认） | 正标量
采样率，指定为正标量。采样率是单位时间内的采样次数。如果时间单位是秒，那么采样率的单位就是赫兹。
# 名称-值对参数
指定可选的以逗号分隔的 Name,Value 参数对。Name 是参数名称，Value 是相应的值。名称必须在引号内。您可以按任意顺序指定多个名称和值对参数，如 Name1,Value1,...,NameN,ValueN。
示例： ty_spectrogram(x; stype ="power", freqrange="centered",OutputTimeDimension="downrows",plotfig=true,freqloc="xaxis")
freqrange - PSD 估计的频率范围"onesided" | "twosided" | "centered"
PSD 估计的频率范围，指定为 "onesided"、"twosided" 或 "centered"。对于实值信号，默认为 "onesided"。对于复值信号，默认为 "twosided"，指定 "onesided" 会导致错误。
"onesided" - 返回实数输入信号的单边频谱图。如果 nfft 为偶数，则 ps 有 nfft/2 + 1 行，并在 [0, π] rad/采样间隔内计算。如果 nfft 为奇数，则 ps 有 (nfft + 1)/2 行，区间为 [0, π] rad/sample。如果指定 fs，那么间隔分别为 [0, fs/2] 周期/单位时间和 [0, fs/2) 周期/单位时间；
"twosided" - 返回真实或复杂信号的双侧频谱图。ps 有 nfft 行，计算区间为 [0, 2π) rad/采样点。如果指定 fs，则区间为 [0, fs) 周期/单位时间；
"centered" - 返回实数或复数信号的居中双面频谱图。如果 nfft 为偶数，则 ps 在 (-π, π] rad/采样间隔内计算。如果 nfft 为奇数，则 ps 在 (-π, π) rad/sample 上计算。如果指定 fs，则间隔分别为 (-fs/2, fs/2] 周期/单位时间和 (-fs/2, fs/2) 周期/单位时间。
stype - 功率谱缩放"psd"（默认） | "power"
功率谱缩放，指定为 "psd "或 "power"。
省略频谱类型或指定 "psd"，将返回功率谱密度；
指定 "power" 时，会根据窗的等效噪声带宽对 PSD 的每个估计值进行缩放。结果就是每个频率的功率估计值。如果使用了 reassign = true 选项，函数会在重新分配之前对每个频率分区的宽度进行 PSD 积分。
reassign - 频谱重分配选项false (默认) | true
功率谱的重分配选项，指定为 true 或 false。 当指定为 true 时对功率谱进行重分配。
数据类型： Bool
MinThreshold - 阈值-Inf（默认） | 实标量
阈值，指定为由 MinThreshold 和以分贝表示的实数标量组成。谱图将 s 的那些元素设置为零，使得 10log10（s）≤thresh。
OutputTimeDimension - 输出时间维度"acrosscolumns"（默认） | "downrows"
输出时间维度，指定为 "acrosscolumns" 或 "downrows" 。如果希望 s 的时间维度沿行，频率维度沿列暑促胡，则将此值设为 "downrows"；如果希望 s 的时间维度按列，频率维度沿行，则将此值设为 "acrosscolumns"。
plotfig - 绘图选项true | false
绘图选项，指定为 true 或 false。 当指定为 true 时绘制图像，否则不绘制。
数据类型： Bool
freqloc - 频率显示轴"yaxis"（默认） | "xaxis"
频率显示轴，指定为 "xaxis" 或 "yaxis"。
"xaxis" - 在 x 轴上显示频率，在 y 轴上显示时间。
"yaxis" -在 y 轴上显示频率，在 x 轴上显示时间。
如果调用 ty_spectrogram输出参数时，则此参数将被忽略。
# 输出参数
s - 短时傅立叶变换矩阵
短时傅里叶变换，以矩阵形式返回。时间从 s 的列开始增加，频率从 0 开始沿行增加。
如果 x 是长度为 Nx 的信号，那么 s 有 k 列，其中：
如果窗是标量，则 k = ⌊（Nx - noverlap）/（window - noverlap）⌋；
如果窗是向量，则 k = ⌊（Nx - noverlap）/（length(window) - noverlap）⌋。
如果 x 为实数且 nfft 为偶数，则 s 有 (nfft/2 + 1) 行；
如果 x 为实数，nfft 为奇数，那么 s 有 (nfft + 1)/2 行；
如果 x 为复数，则 s 有 nfft 行。
s 不受 reassign 选项的影响。
w - 归一化频率向量
归一化频率，以向量形式返回。w 的长度等于 s 的行数。
t - 时间瞬时向量
时间瞬时，以向量形式返回。t 中的时间值对应于每个片段的中点。
f - 周期频率向量
周期频率，以向量形式返回。f 的长度等于 s 的行数。
ps - 功率谱密度或功率谱矩阵
功率谱密度（PSD）或功率谱，以矩阵形式返回。
如果 x 为实数，则 ps 包含每个片段的 PSD 或功率谱的单边修正周期图估计值；
如果 x 为复数，或者您指定了一个频率向量，那么 ps 将包含每个分段的 PSD 或功率谱的双面修正周期图估计值。
fc, tc - 能量中心频率和时间矩阵
能量中心频率和时间，以与短时傅里叶变换相同大小的矩阵形式返回。如果没有指定采样率，则 fc 的元素将以归一化频率的形式返回。
fc 和 tc 仅在输入 reassign = true 时返回。
# 参考文献
[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.
[2] Rabiner, Lawrence R., and Ronald W. Schafer. Digital Processing of Speech Signals. Englewood Cliffs, NJ: Prentice-Hall, 1978.
[3] Chassande-Motin, Éric, François Auger, and Patrick Flandrin. “Reassignment.” In Time-Frequency Analysis: Concepts and Methods. Edited by Franz Hlawatsch and François Auger. London: ISTE/John Wiley and Sons, 2008.
[4] Fulop, Sean A., and Kelly Fitz. “Algorithms for computing the time-corrected instantaneous frequency (reassigned) ty_spectrogram, with applications.” Journal of the Acoustical Society of America. Vol. 119, January 2006, pp. 360–371.