从相关函数到相干函数：语音信号的时频相关分析入门

在语音、回声消除、反馈抑制、阵列处理这些问题里，我们经常会说两个信号“相关”“相干”“谱上有关联”。这些词听起来很接近，但真正落到分析时，它们对应的是不同指标。本文用一段类语音信号作为例子，讲解**自相关函数、自功率谱、互相关函数、互功率谱、幅值平方相干函数 ** 之间的联系与区别。

本文同步发在个人博客网站VoxWorking · 声学与音频知识平台，网站定期更新音频处理、语音处理和声学相关知识，欢迎大家浏览交流~

1. 实验信号

图中上面是原始类语音信号 x(n)，下面是延迟和滤波后的观测信号 y(n)。可以看到两者有相似的起伏结构，但发生了时间偏移、频谱改变和噪声污染。

2. 自相关函数：信号和自己的相似程度

自相关函数描述一个信号和它自身延迟版本之间的相似程度：

$$ R_{xx}[k] = E[x(n)x(n-k)] $$

其中 k 是 lag，表示延迟多少个采样点。如果 Rxx[k] 在 k=0 以外仍然很大，说明信号不是白噪声，而是存在可预测结构。语音通常会有明显自相关，因为它包含短时谱包络相关、基音和谐波导致的周期相关、包络变化导致的低频调制等。在 MATLAB 里可以这样计算，其中 'coeff' 表示归一化，结果会落在大约 [-1, 1] 的范围内：

maxLag = round(0.08 * Fs);          % 分析 +/-80 ms
[rxx, lags] = xcorr(x, maxLag, 'coeff');
plot(lags/Fs*1000, rxx);
xlabel('Lag (ms)');
ylabel('Normalized autocorrelation');
grid on;

3. 自功率谱：自相关的频域版本

自功率谱密度描述信号能量在频率上的分布，常写成：

$$ S_{xx}(\omega) $$

它和自相关函数是一对傅里叶变换：

$$ S_{xx}(\omega) = \sum_{k=-\infty}^{\infty} R_{xx}[k]e^{-j\omega k} $$

反过来也有：

$$ R_{xx}[k] = \frac{1}{2\pi}\int_{-\pi}^{\pi}S_{xx}(\omega)e^{j\omega k}d\omega $$

这就是常说的 Wiener-Khinchin 关系。自功率谱和自相关有以下关系：（1）频谱平坦，时域上更接近不相关；（2）频谱起伏明显，时域上通常存在相关。（3）谐波峰明显，时域上容易出现周期相关峰。

上图是 x(n) 的自相关，下图是它的功率谱。可以看到语音不是白噪声：自相关在非零 lag 上仍然明显，功率谱也不是平坦的，而是有强弱不均的谱结构。MATLAB 中可以用 pwelch 估计功率谱：

nfft = 2048;
win = hann(1024, 'periodic');
noverlap = 512;

[Pxx, f] = pwelch(x, win, noverlap, nfft, Fs);
semilogy(f, Pxx);
xlim([0 5000]);
grid on;

4. 互相关函数：两个信号在哪个延时上最像

互相关函数描述两个信号之间的延迟相似性。这里可以注意下互相关 lag 正负号的定义。为了让图更符合工程直觉，本文采用下面这个约定：正 lag 表示 y(n) 相对 x(n) 滞后。也就是说，如果： $$ y(n)=x(n-D) $$

那么互相关峰值应该出现在：

$$ k \approx +D $$

对应的互相关可以理解为：

$$ R_{yx}[k] = E[y(n+k)x(n)] $$

在 MATLAB 里，为了得到这个 lag 方向，代码使用：

[ryx, lags] = xcorr(y, x, maxLag, 'coeff');
plot(lags/Fs*1000, ryx);

如果改成 xcorr(x, y)，同一个物理延时会出现在相反方向。这不是信号变了，而是互相关定义的正负号变了。

5. 互功率谱：互相关的频域版本

互功率谱，也叫 cross power spectral density，常写作：

$$ S_{yx}(\omega) $$

它和互相关函数同样是一对傅里叶变换：

$$ S_{yx}(\omega)=\sum_{k=-\infty}^{\infty}R_{yx}[k]e^{-j\omega k} $$

互功率谱是复数，包含两部分信息：（1）幅度：两路信号在某个频率上共同变化有多强。（2）相位：两路信号在该频率上的相位差。如果 y(n) 近似是 x(n) 延迟 D 点得到：

$$ y(n) \approx x(n-D) $$

那么频域中近似有：

$$ Y(\omega)\approx X(\omega)e^{-j\omega D} $$

因此互功率谱相位会出现近似线性斜率。注意，相位正负号同样取决于你用的是 Sxy 还是 Syx，所以工程上更重要的是保持互相关、互谱、代码和图注的顺序一致：

$$ \angle S_{yx}(\omega)\approx +\omega D $$

这也是为什么延时估计不仅可以看互相关峰，也可以看互谱相位斜率。MATLAB 中可以用 cpsd 估计互功率谱：

[Pyx, f] = cpsd(y, x, win, noverlap, nfft, Fs);

subplot(2,1,1);
semilogy(f, abs(Pyx));
xlim([0 5000]);
grid on;

subplot(2,1,2);
plot(f, unwrap(angle(Pyx)));
xlim([0 3000]);
grid on;
上图中，互相关在约 +35 ms 附近出现明显峰值，说明 `y(n)` 中确实含有 `x(n)` 的延迟版本。下图是互功率谱幅度，说明两路信号在哪些频率上存在共同成分。

6. 幅值平方相干函数：归一化后的频域相关性

幅值平方相干函数，常写作 magnitude-squared coherence，定义为：

$$ \gamma_{xy}^2(\omega)= \frac{|S_{xy}(\omega)|^2}{S_{xx}(\omega)S_{yy}(\omega)} $$

它的取值范围通常是：

$$ 0 \le \gamma_{xy}^2(\omega) \le 1 $$

可以这样理解：

注意，相干函数不是互功率谱。它是归一化后的指标，互功率谱 Sxy 或 Syx 会受能量大小影响；相干函数把 Sxx 和 Syy 的能量也归一化了。

上图是互功率谱相位。由于 y(n) 是 x(n) 延迟和滤波后得到的，所以低频到中频范围内能看到较明显的相位趋势。相位曲线向上还是向下，取决于你计算的是 Sxy 还是 Syx；本文代码统一使用 Syx，与互相关正 lag 表示延迟的约定保持一致。下图是相干函数，可以看到某些频带相干性较强，而噪声较多或能量较弱的频带相干性下降。MATLAB 中可以用 mscohere：

[Cxy, fcoh] = mscohere(x, y, win, noverlap, nfft, Fs);
plot(fcoh, Cxy);
ylim([0 1]);
xlim([0 5000]);
grid on;

7. 互谱变了，相干不一定变

如果两路信号同时经过同一个线性时不变滤波器 A(z)：

$$ x_w(n)=A(z)x(n) $$

$$ y_w(n)=A(z)y(n) $$

那么：

$$ S_{x_wy_w}(\omega)=|A(\omega)|^2S_{xy}(\omega) $$

所以互功率谱的幅度会被改变。但相干函数变成：

$$ \gamma_{x_wy_w}^2(\omega)= \frac{||A(\omega)|^2S_{xy}(\omega)|^2} {|A(\omega)|^2S_{xx}(\omega)\cdot |A(\omega)|^2S_{yy}(\omega)} $$

化简后：

$$ \gamma_{x_wy_w}^2(\omega)= \gamma_{xy}^2(\omega) $$

也就是说，在理想条件下，两路经过同一个 LTI 滤波器，相干函数可以保持不变。我们平时说“相关性下降”，可能指的是：互相关峰值下降、互谱幅度下降，但它不一定等价于相干函数下降。

8. 小结

自相关和自功率谱是一对时频描述，描述的是单个信号自身的结构。互相关和互功率谱也是一对时频描述，描述的是两个信号之间的延时关系和频域共同成分。幅值平方相干函数则是在互功率谱基础上做归一化，用来衡量某个频率上两路信号的线性相关程度。

对语音信号来说，这些指标很有用，因为语音同时具有短时谱包络、基音周期、谐波结构和非平稳包络。只看波形往往很难判断问题来源，而把相关函数、功率谱、互谱和相干函数放在一起看，很多隐藏结构就会变得清楚。

后面文章我们会讨论：语音中的短时相关、周期相关、长时相关分别长什么样，以及它们为什么会影响闭环系统中的自适应估计。

9. 实验代码

clear; close all; clc;
rng(7);

scriptDir = fileparts(mfilename('fullpath'));
if isempty(scriptDir)
    scriptDir = pwd;
end
figDir = fullfile(scriptDir, '..', 'figures');
if ~exist(figDir, 'dir')
    mkdir(figDir);
end

Fs = 16000;
dur = 2.5;
t = (0:round(Fs * dur)-1)' / Fs;
lineWidth = 1.6;
titleFontSize = 17;
labelFontSize = 15;
tickFontSize = 13;

%% Build a simple voice-like signal
f0 = 165 + 18 * sin(2*pi*0.55*t) + 8 * sin(2*pi*1.2*t);
phase = 2*pi*cumsum(f0) / Fs;

x = zeros(size(t));
for h = 1:22
    x = x + (1/h) * sin(h * phase + 0.13*h);
end

env = (0.55 + 0.35 * sin(2*pi*1.1*t).^2) .* ...
      (0.8 + 0.2 * sin(2*pi*0.25*t + 0.6));
x = x .* env;

x = resonator_filter(x, Fs, 700, 0.965);
x = resonator_filter(x, Fs, 1250, 0.955);
x = resonator_filter(x, Fs, 2600, 0.94);
x = x + 0.015 * randn(size(x));
x = x / max(abs(x) + eps);

%% Build an observed signal y(n): delay + filtering + noise
delayMs = 35;
D = round(delayMs * 1e-3 * Fs);
y = zeros(size(x));
y(D+1:end) = 0.82 * x(1:end-D);
y = resonator_filter(y, Fs, 1600, 0.92);
y = y + 0.05 * randn(size(y));
y = y / max(abs(y) + eps);

%% Correlation analysis
maxLagMs = 80;
maxLag = round(maxLagMs * 1e-3 * Fs);

[rxx, lagXx] = xcorr(x, maxLag, 'coeff');

% Delay-estimation convention used in this article:
% positive lag means y is delayed relative to x.
% Therefore, for y(n) = x(n-D), xcorr(y, x) peaks near +D.
[ryx, lagYx] = xcorr(y, x, maxLag, 'coeff');

%% Spectral analysis
nfft = 2048;
win = hann(1024, 'periodic');
noverlap = 512;

[Pxx, f] = pwelch(x, win, noverlap, nfft, Fs);
[Pyy, ~] = pwelch(y, win, noverlap, nfft, Fs);
[Pyx, fc] = cpsd(y, x, win, noverlap, nfft, Fs);
[Cxy, fcoh] = mscohere(x, y, win, noverlap, nfft, Fs);

%% Plot results
figure('Name', 'Signals');
subplot(2,1,1);
plot(t(1:round(0.22*Fs))*1000, x(1:round(0.22*Fs)), 'LineWidth', lineWidth);
grid on; set(gca, 'FontSize', tickFontSize);
title('voice-like signal x(n)', 'FontSize', titleFontSize);
ylabel('Amplitude', 'FontSize', labelFontSize);
subplot(2,1,2);
plot(t(1:round(0.22*Fs))*1000, y(1:round(0.22*Fs)), 'LineWidth', lineWidth);
grid on; set(gca, 'FontSize', tickFontSize);
title(sprintf('observed signal y(n): delayed by %d ms', delayMs), 'FontSize', titleFontSize);
xlabel('Time (ms)', 'FontSize', labelFontSize); ylabel('Amplitude', 'FontSize', labelFontSize);
saveas(gcf, fullfile(figDir, 'fig1_signals_time_domain.png'));

figure('Name', 'Autocorrelation and PSD');
subplot(2,1,1);
plot(lagXx/Fs*1000, rxx, 'LineWidth', lineWidth); grid on;
set(gca, 'FontSize', tickFontSize);
title('Autocorrelation R_{xx}[k]', 'FontSize', titleFontSize);
xlabel('Lag (ms)', 'FontSize', labelFontSize); ylabel('Normalized correlation', 'FontSize', labelFontSize);
subplot(2,1,2);
plot(f, 20*log10(Pxx + eps), 'LineWidth', lineWidth); grid on; xlim([0 5000]);
set(gca, 'FontSize', tickFontSize);
title('Power spectral density S_{xx}(f)', 'FontSize', titleFontSize);
xlabel('Frequency (Hz)', 'FontSize', labelFontSize); ylabel('Magnitude dB', 'FontSize', labelFontSize);
saveas(gcf, fullfile(figDir, 'fig2_autocorr_psd.png'));

figure('Name', 'Cross-correlation and CPSD');
subplot(2,1,1);
plot(lagYx/Fs*1000, ryx, 'LineWidth', lineWidth); grid on; hold on;
delayLine = xline(delayMs, '--r', 'expected delay', 'LineWidth', lineWidth);
if isprop(delayLine, 'FontSize')
    delayLine.FontSize = tickFontSize;
end
set(gca, 'FontSize', tickFontSize);
title('Cross-correlation R_{yx}[k]', 'FontSize', titleFontSize);
xlabel('Lag (ms)', 'FontSize', labelFontSize); ylabel('Normalized correlation', 'FontSize', labelFontSize);
subplot(2,1,2);
plot(fc, 20*log10(abs(Pyx) + eps), 'LineWidth', lineWidth); grid on; xlim([0 5000]);
set(gca, 'FontSize', tickFontSize);
title('|Cross power spectral density S_{yx}(f)|', 'FontSize', titleFontSize);
xlabel('Frequency (Hz)', 'FontSize', labelFontSize); ylabel('Magnitude dB', 'FontSize', labelFontSize);
saveas(gcf, fullfile(figDir, 'fig3_crosscorr_cpsd.png'));

figure('Name', 'CPSD phase and coherence');
subplot(2,1,1);
plot(fc, unwrap(angle(Pyx)), 'LineWidth', lineWidth); grid on; xlim([0 3000]);
set(gca, 'FontSize', tickFontSize);
title('CPSD phase, angle(S_{yx})', 'FontSize', titleFontSize);
ylabel('Unwrapped phase (rad)', 'FontSize', labelFontSize);
subplot(2,1,2);
plot(fcoh, Cxy, 'LineWidth', lineWidth); grid on; xlim([0 5000]); ylim([0 1.05]);
set(gca, 'FontSize', tickFontSize);
title('Magnitude-squared coherence', 'FontSize', titleFontSize);
xlabel('Frequency (Hz)', 'FontSize', labelFontSize); ylabel('Coherence', 'FontSize', labelFontSize);
saveas(gcf, fullfile(figDir, 'fig4_cpsd_phase_coherence.png'));

%% Local helper
function y = resonator_filter(x, Fs, f0, radius)
    theta = 2*pi*f0/Fs;
    a1 = 2 * radius * cos(theta);
    a2 = -(radius^2);
    y = zeros(size(x));
    for n = 1:length(x)
        y(n) = x(n);
        if n >= 2
            y(n) = y(n) + a1 * y(n-1);
        end
        if n >= 3
            y(n) = y(n) + a2 * y(n-2);
        end
    end
end