# ty_mscohere
幅度平方相干性
函数库: TySignalProcessing
# 语法
cxy, = ty_mscohere(x, y)
cxy, = ty_mscohere(x, y, window)
cxy, = ty_mscohere(x, y, window, noverlap)
cxy, = ty_mscohere(x, y, window, noverlap, nfft)
cxy, = ty_mscohere(___; mimo = true)
cxy, w = ty_mscohere(___)
cxy, f = ty_mscohere(___, fs)
cxy, w = ty_mscohere(x, y, window, noverlap, w)
cxy, f = ty_mscohere(x, y, window, noverlap, f, fs)
___ = ty_mscohere(x, y, ___; freqrange = String)
ty_mscohere(___; plotfig = true)
# 说明
cxy, = ty_mscohere(x, y) 找到输入信号 x 和 y 的幅度平方相干估计值 cxy。
如果 x 和 y 都是向量，它们必须有相同的长度；
如果其中一个信号是矩阵，另一个是向量，那么向量的长度必须等于矩阵的行数。该函数对向量进行扩展，并返回一个逐列幅度平方相干估计的矩阵；
如果 x 和 y 是行数相同但列数不同的矩阵，那么 ty_mscohere 会返回一个多重相干矩阵。cxy 的第 m 列包含所有输入信号和第 m 个输出信号之间的相关程度的估计；
如果 x 和 y 是大小相等的矩阵，那么 ty_mscohere 是逐列操作的：cxy[:, n], = ty_mscohere(x[:, n], y[:, n])。要获得一个多重相干矩阵，请在参数中输入 mimo = true。
cxy, = ty_mscohere(x, y, window) 使用 window 对 x 和 y 分段，并进行加窗处理。至少需要有两个分段。否则，所有频率的幅度-平方相干性为 1。在 MIMO 情况下，段的数量必须大于输入通道的数量。
cxy, = ty_mscohere(x, y, window, noverlap) 使用相邻段之间重叠的 noverlap 样本。
cxy, = ty_mscohere(x, y, window, noverlap, nfft) 使用 nfft 采样点来计算离散傅里叶变换。
cxy, = ty_mscohere(___; mimo = true) 计算矩阵输入的个多重相干矩阵。这种语法可以包含此前语法中输入参数的任何组合。
cxy, w = ty_mscohere(___) 返回一个归一化的频率向量 w，在这个频率上的幅度平方相干性估计。
cxy, f = ty_mscohere(___, fs) 返回一个频率的向量 f，用采样率 fs 表示，在这个采样率上估计幅度平方相干性。fs 必须是 ty_mscohere 的第六个数字输入。要输入一个采样率，并仍然使用前面的可选参数的默认值，请将这些参数指定为空 [ ]。
cxy, w = ty_mscohere(x, y, window, noverlap, w) 返回 w 中指定的归一化频率下的幅度平方相干性估计。
cxy, f = ty_mscohere(x, y, window, noverlap, f, fs) 返回 f 中指定的频率下的幅度平方相干性估计。
___ = ty_mscohere(x, y, ___; freqrange = String) 返回由 freqrange 指定的频率范围内的幅度平方相干性估计。freqrange 的有效选项是 "onesided"、"twosided" 和 "centered"。
ty_mscohere(___; plotfig = true) 在没有输出参数的情况下，在当前图形窗口中绘制出幅度平方相干的估计值。
# 示例
多重相干与普通相干
产生一个随机的双通道信号 x，通过低通滤波产生另一个信号 y，并将它们加在一起。指定一个截止频率为 0.3 π 的 30 阶 FIR 滤波器，使用矩形窗口设计。
using TyPlot
using TySignalProcessing
using TyMath
rng = MT19937ar(1234)
h = fir1(30, 0.3, rectwin(31))
x = randn(rng, 16384, 2)
y = sum(filter1(h, 1, x)[1]; dims=2)
计算 x 和 y 的多重相干估计值。用 1024 个样本的 Hann 窗口对信号进行处理。指定 512 个相邻段之间的重叠样本和 1024 个 DFT 点。绘制估计值。
noverlap = 512
nfft = 1024
ty_mscohere(x, y, hann(nfft), noverlap, nfft; mimo=true, plotfig=true)
将相干性估计值与滤波器的频率响应进行比较。相干性的下降对应于频率响应的零点。
H, f = freqz(h)
hold("on")
yyaxis("right")
ylim([-80, 10])
plot(f ./ pi, 20 * log10.(abs.(H)))
hold("off")
计算并绘制 x 和 y 的幅度平方相干估计值。
figure()
ty_mscohere(x, y, hann(nfft), noverlap, nfft; plotfig=true)
MIMO 系统的相干性
产生两个多通道信号，每个信号的采样频率为 1 kHz，持续 2 秒。第一个信号，即输入信号，由三个频率为 120 Hz、360 Hz 和 480 Hz 的正弦波组成。第二个信号，即输出，由两个频率为 120 Hz 和 360 Hz 的正弦波组成。其中一个正弦波比第一个信号滞后 π/2。另一个正弦信号的滞后为 π/4。这两个信号都被嵌入到白高斯噪声中。
using TyPlot
using TySignalProcessing
using TyMath
rng = MT19937ar(1234)
fs = 1000
f = 120
t = (0:(1 / fs):(2 - 1 / fs))
inpt = sin.(2 * pi * f * [1 3 4] .* t)
inpt = inpt + randn(rng, size(inpt))
oupt = sin.(2 * pi * f * [1 3] .* t .- [pi / 2 pi / 4])
oupt = oupt + randn(rng, size(oupt))
估计所有输入信号和每个输出通道之间的相关程度。ty_mscohere 为每个输出通道返回一个相干函数。相干函数在输入和输出共享的频率上达到最大值。
h = hamming(100)
Cxy, f = ty_mscohere(inpt, oupt, h, [], [], fs; mimo=true)
for k in axes(oupt, 2)
subplot(size(oupt, 2), 1, k)
plot(f, Cxy[:, k])
title("Output" * string(k) * ", All Inputs")
end
tightlayout()
切换输入和输出信号，计算多重相干函数。使用相同的汉明窗。在 480 Hz 时，输入和输出之间没有关联性。因此，在第三个相关函数中没有峰值。
Cxy, f = ty_mscohere(oupt, inpt, hamming(100), [], [], fs; mimo=true)
figure()
for k in axes(inpt, 2)
subplot(size(inpt, 2), 1, k)
plot(f, Cxy[:, k])
title("Input " * string(k) * ", All Outputs")
end
tightlayout()
重复计算，使用 ty_mscohere 的绘图功能。
figure()
ty_mscohere(oupt, inpt, hamming(100), [], [], fs; mimo=true, plotfig=true)
计算第二信号和第一信号的前两个通道的普通相干函数。非峰值部分与多重相干函数不同。
Cxy, f = ty_mscohere(oupt, inpt[:, 1:2], hamming(100), [], [], fs)
figure()
plot(f, Cxy)
通过计算最大相干点的互谱的角度来找到相位差。
Pxy, = ty_cpsd(oupt, inpt[:, 1:2], hamming(100), [], [], fs)
mxx = argmax(abs.(Cxy); dims=1)
for k in 1:2
@printf("Phase lag %d = %5.2f*pi\n", k, angle(Pxy[mxx[k]]) / pi)
end
Phase lag 1 = -0.52*pi
Phase lag 2 = -0.27*pi
修改幅度平方相干图
产生两个正弦信号，每个信号以 1 kHz 采样 1 秒。每个正弦信号的频率为 250 Hz。其中一个信号在相位上比另一个滞后 π/3 弧度。将两个信号嵌入单位方差的高斯白噪声中。
using TyPlot
using TySignalProcessing
using TyMath
rng = MT19937ar(1234)
fs = 1000
f = 250
t = 0:(1 / fs):(1 - 1 / fs)
um = sin.(2 * pi * f * t) + rand(rng, size(t))
un = sin.(2 * pi * f * t .- pi / 3) + rand(rng, size(t))
使用 ty_mscohere 来计算和绘制信号的幅度-平方相干性。
ty_mscohere(um, un, [], [], [], fs; plotfig=true)
修改绘图的标题，x 轴的标签，以及 y 轴的界限。
title("Magnitude-Squared Coherence")
xlabel("f (Hz)")
ylim([0 1.1])
使用 gca 获取当前坐标轴。改变刻度线的位置。删除 y 轴的标签。
ax = gca()
xticks(0:250:500)
yticks(0:0.25:1)
ylabel("")
调用 get_lines 属性来改变绘图线的颜色和宽度。
line = ax.get_lines()
line[1].set_color([0.8 0 0])
line[1].set_linewidth(1.5)
或者，使用其他方式更改属性。
line[1].set_color([0 0.4 0])
line[1].set_linestyle("--")
line[1].set_linewidth(1)
# 输入参数
x, y - 输入信号向量 | 矩阵
输入信号，指定为向量或矩阵。
数据类型： Float | Int
window - 窗整数 | 向量 | [ ]
窗，指定为一个整数或一个行或列向量。使用 window 将信号分成若干段：
如果 window 是一个整数，那么 ty_mscohere 将 x 和 y 分成长度为 window 的片段，并用该长度的汉明窗对每个片段进行处理；
如果 window 是一个向量，那么 ty_mscohere 将 x 和 y 分成与向量相同长度的片段，并使用 window 对每个片段进行窗处理。
如果 x 和 y 的长度不能精确地划分为整数个段，并有 noverlap 重叠的样本，那么信号就会被相应地截断。
如果你指定 window 为 [ ]，那么 ty_mscohere 使用一个 Hamming 窗，这样 x 和 y 就被分成 8 段，并有 noverlap 重叠的样本。
数据类型： Int | Float
noverlap - 重叠的样本数正整数 | [ ]
DFT 点的数量，指定为一个正整数。如果你指定 nfft 为空，那么 ty_mscohere 将这个参数设置为 maximum(256, 2p)，其中 p = ⌈log2 N⌉，用于长度为 N 的输入信号。
数据类型： Int
nfft - DFT 点的数量正整数 | [ ]
DFT 点的数量，指定为一个正整数。如果你指定 nfft 为空，那么 ty_mscohere 将这个参数设置为
数据类型： Int
fs - 采样率正标量
采样率，指定为一个正标量。采样率是每单位时间的采样数。如果时间的单位是秒，那么采样率的单位是 Hz。
数据类型： Float | Int
w - 归一化的频率向量
归一化的频率，指定为至少有两个元素的行或列向量。归一化频率的单位是 rad/sample。
数据类型： Float | Int
f - 频率向量
频率，指定为一个至少有两个元素的行或列向量。频率的单位是每单位时间的周期。单位时间是由采样率 fs 指定的。如果 fs 的单位是 samples/second，那么 f 的单位就是 Hz。
数据类型： Float | Int
# 名称-值参数
指定可选的、以逗号分隔的 Key=Value 对组参数。Key 为参数名称，Value 为对应的值。您可采用任意顺序指定多个名称-值对组参数，如 Key1=Value1,...,KeyN=ValueN 所示。
示例: ty_mscohere(x,y; freqrange="twosided", mimo=true, plotfig=true)
freqrange - 幅度平方相干估计的频率范围"onesided" | "twosided" | "centered"
幅度平方相干估计的频率范围，指定为 "onesided"、"twosided" 和 "centered" 中的一个。对于实值信号，默认为 "onesided"，对于复值信号，默认为 "twosided"。
"onesided" - 返回两个实值输入信号 x 和 y 之间的单边相干估计值；
- 如果 nfft 是偶数，cxy 有 (nfft + 1) / 2 行，区间是 [0, π] rad/sample；
- 如果 nfft 是奇数，cxy 有 (nfft + 1) / 2 行，区间为 [0, π) rad/sample；
- 如果你指定 fs，对于偶数 nfft，相应的间隔是 [0, fs/2] cycles/unit 时间，对于奇数 nfft，相应的间隔是 [0, fs/2) cycles/unit 时间。
"twosided" - 返回两个实值或复值输入信号 x 和 y 之间的幅度平方相干估计值的双边估计。在这种情况下，cxy 有 nfft 行，并在区间 [0, 2π) rad/sample 上计算。如果你指定 fs，那么这个区间就是 [0, fs] cycles/unit 时间；
"centered" - 返回两个实值或复值输入信号 x 和 y 之间的幅度平方相干估计的中心双边估计值。如果你指定了 fs，那么对于偶数 nfft 来说，相应的区间是 (-fs/2, fs/2] cycles/unit 时间，对于奇数 nfft 来说，是 (-fs/2, fs/2) cycles/unit 时间。
数据类型： String
mimo - 多输入/多输出选项true | false
多输入/多输出选项，指定为 true 或 false。 当指定为 true 时计算多输入多输出的幅度平方相干性。
数据类型： Bool
plotfig - 绘图选项true | false
绘图选项，指定为 true 或 false。 当指定为 true 时绘制图像，否则不绘制。
数据类型： Bool
# 输出参数
cxy - 幅度平方相干估计值向量 | 矩阵
幅度平方相干估计值，以向量或矩阵形式返回。
数据类型： Float
w - 归一化的频率向量
归一化的频率，作为实值列向量返回。
数据类型： Float | Int
f - 频率向量
频率，作为一个实值列向量返回。
数据类型： Float | Int
# 更多关于
幅度平方相干性
幅度平方相干性估计是一个频率的函数，其值在 0 和 1 之间。这些值表示 x 在每个频率下与 y 的对应程度。幅度平方相干性是 x 和 y 的功率谱密度
对于多输入/多输出系统来说，多重相干函数变为:
为第 i 个输出信号，其中：
X 对应的是 m 个输入的数组;
是输入和 yi 之间的交叉功率谱密度的 m 维向量; 是输入的功率谱密度和交叉功率谱密度的 m×m 矩阵; 是输出的功率谱密度; dagger(†) 代表复数共轭转置。
# 算法
ty_mscohere 使用 Welch 的重叠平均周期图方法 [3]，[5] 来估计幅度-平方相干函数 [2]。
# 参考文献
[1] Gómez González, A., J. Rodríguez, X. Sagartzazu, A. Schumacher, and I. Isasa. “Multiple Coherence Method in Time Domain for the Analysis of the Transmission Paths of Noise and Vibrations with Non-Stationary Signals.” Proceedings of the 2010 International Conference of Noise and Vibration Engineering, ISMA2010-USD2010. pp. 3927–3941.
[2] Kay, Steven M. Modern Spectral Estimation. Englewood Cliffs, NJ: Prentice-Hall, 1988.
[3] Rabiner, Lawrence R., and Bernard Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975.
[4] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.
[5] Welch, Peter D. “The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms.” IEEE® Transactions on Audio and Electroacoustics. Vol. AU-15, 1967, pp. 70–73.