简介:本资源是一套面向信号处理初学者与MATLAB实践者的教学辅助代码包,聚焦自相关、协方差、信号方差等核心统计特性分析,解决理论理解抽象、公式推导与实际计算脱节的学习痛点,适用于通信、电子、自动化等专业本科生及工程入门者。压缩包为RAR格式,共含4个MATLAB脚本文件(.m),总大小仅2KB,轻量精炼,涵盖自相关函数实现、信号方差计算、自协方差与互协方差仿真等关键模块,每个脚本均对应典型信号处理场景(如平稳性检验、延迟相关性分析、噪声特性评估),便于逐行调试与结果可视化。已有388人学习下载,资源虽小但结构完整:代码注释清晰、变量命名规范、含典型测试信号(如正弦加噪序列)与绘图指令,可直接运行观察时域相关性衰减趋势、协方差矩阵对称性等本质特征,是夯实信号统计分析基础的实用工具脚本集。 刚拿到一个Autocorrelation_function.rar的压缩包时,我第一反应是:这又是个随手存的信号处理作业。但解压打开后,里面的文件名把"信号方差、协方差、相关函数、自相关"这些关键词全串在了一起,我突然意识到这其实不是零散代码,而是一条完整的信号统计分析链路。在 MATLAB 里搞信号处理的同学,绝大多数都用过xcorr、cov、var这些函数,但真要说清楚"自相关和协方差到底什么关系""为什么功率信号的自相关要归一化""互相关的峰值为什么能用来算时延",不少人会卡壳。这篇就把这条链路完整捋一遍,从概念到 MATLAB 实现,再到实际场景中的参数取舍和踩坑记录,给准备用相关分析处理信号的你一份能直接照着用的参考资料。无论你是刚接触随机信号分析的学生,还是需要做周期检测、时延估计的工程师,都能从中找到对应自己问题的答案。
1. 自相关函数到底在算什么:从一个被误读最多的概念说起
先说一个我在论坛上见过无数次的误解:很多人把自相关当成"信号和自身的相似度",这句话对了一半,但它忽略了一个关键的数学事实——自相关本质上是一个统计量,不是单纯逐点比较的相似度。它衡量的是"信号与其自身延迟 k 个采样点之后,这两个序列之间的线性相关程度"。这在随机信号分析里意义重大,因为对于平稳随机过程,自相关函数直接反映了信号在不同时间偏移下的统计依赖特性。
1.1 从协方差到自相关的两条推导路径
要理解自相关,可以先从协方差入手。协方差描述的是两个随机变量 X 和 Y 一起变化的程度,公式是:
cov(X, Y) = E[(X - μX)(Y - μY)]如果令 Y 等于 X 自身向后延迟 k 个时刻的序列,也就是Y(t) = X(t + k),那么协方差就变成了"信号与自身延迟序列的协方差",这时它被命名为自协方差函数:
γ(k) = E[(X(t) - μ)(X(t + k) - μ)]而自相关函数则有两种常见的定义口径。第一种是直接把自协方差归一化,得到的是数值范围在 [-1, 1] 的相关系数序列:
ρ(k) = γ(k) / γ(0)第二种是信号处理领域更常用的定义,不对均值做中心化处理,直接计算两个序列乘积的期望:
R(k) = E[X(t) · X(t + k)]这两种定义各有适用场景。前者在统计学和随机过程理论中更常见,用于分析信号内部的线性依赖结构;后者在工程实现中更常见,比如功率谱估计里,维纳-辛钦定理说的就是功率谱密度与自相关函数互为傅里叶变换对,这里的自相关用的就是非中心化的R(k)。
1.2 为什么要把"信号方差"和自相关绑在一起看
标题里出现了"信号方差",这个组合其实暗含了信号分析中一个非常重要的关系:自相关函数在零延迟处的值,就等于信号的方差(对于零均值信号),或者等于信号的均方值(对于非零均值信号)。这一点在工程上极其有用。
举个例子,R(0) = E[X(t)²],对于零均值信号来说,E[X(t)²]就是方差。所以当你算完一条自相关曲线,看R(0)这一个点,就能得到信号的功率信息。而R(k)(k 不等于 0)则告诉你信号在不同时间偏移下的统计相关性有多强。对于一个纯随机白噪声,除了R(0)之外,其他延迟处的自相关值都趋近于 0,因为白噪声不同时刻的样本之间毫无关联;对于正弦波这类周期信号,自相关会在延迟等于周期的整数倍处出现峰值,而且这个峰值不会随延迟增大而衰减——这个特性正是用自相关做周期检测的数学基础。
2. MATLAB 里自相关与协方差计算的主干函数:xcorr、xcov、autocorr 怎么选
MATLAB 提供的相关函数其实不止一个,容易让人选择困难。xcorr、xcov、autocorr、corrcoef、cov,五六个函数,名字相近,输出格式和默认行为又各有差别,混用起来很容易翻车。这里把它们的本质区别和使用场景理清楚。
2.1 xcorr 与 xcov:一字之差,差了一个"去均值"
xcorr用于计算互相关或自相关,默认采用非中心化定义。xcov用于计算互协方差或自协方差,它会先减去各自信号的均值再计算相关。也就是说,xcov(x, x)的结果与xcorr(x - mean(x), x - mean(x))是等价的。
对于非零均值信号,这两者的差别是实质性的,不是小修小补。假设你分析一段带有直流偏置的振动传感器数据,如果用xcorr,直流分量会把相关函数整体抬高,导致你在零延迟处看到一个大尖峰,容易掩盖掉信号本身的周期性结构。而xcov去掉了均值,能够更干净地反映出波动部分的相关性。
实际使用时,我的建议是:分析周期性和波动特征优先用xcov或对信号先做去均值处理,再用xcorr;分析信号能量分布或需要配合功率谱计算时,用xcorr。
2.2 autocorr:让结果更接近统计学惯例的封裝
如果安装了 Econometrics Toolbox,可以使用autocorr函数。它返回的是归一化的自相关序列,即ρ(k) = γ(k) / γ(0),范围严格落在 [-1, 1] 内,并且默认会画出带置信区间的棒图。这个函数非常适合做时间序列分析,但在纯信号处理任务里,它的封装程度太高,返回的 lag 比较难直接映射到物理延迟,而且不能像xcorr那样计算两个不同信号之间的互相关。
下面用一个简单例子对比三个函数的输出形态。假设生成一段 1000 点的正弦信号加噪声:
fs = 1000; t = (0:999) / fs; x = sin(2*pi*50*t) + 0.3*randn(1, 1000); % xcorr 非中心化 [r1, lags1] = xcorr(x, 'coeff'); % xcov 中心化 [r2, lags2] = xcov(x, 'coeff'); % autocorr 统计学归一化 [r3, lags3] = autocorr(x, 'NumLags', 100);这三个结果在零延迟处都为 1,但在远离零延迟的位置会有细微差别,尤其是当信号存在直流分量时,差别会非常明显。
2.3 输出数组的长度问题:为什么算完自相关,数据量翻倍了
新手最容易懵的地方在于:输入一段长度 N 的信号,xcorr输出的数组长度是2N - 1。原因是互相关计算要把第二个信号相对于第一个信号进行正负两个方向的平移,所以延迟范围是-(N-1)到N-1。lags数组的中间点(第 N 个元素)对应延迟 0,这一点一定要记牢,否则取峰值位置时很容易取错,我在延迟估计里没少吃这个亏。
如果你想限制延迟范围,xcorr提供了maxlag参数:
[r, lags] = xcorr(x, y, 200, 'coeff');这样只计算延迟从 -200 到 200 的相关值,输出长度为 401,计算量大幅减少,峰值定位也更直观。做时延估计时建议加 maxlag,既能省时间,又能避免找到错误周期对应的远距离峰值。
3. 归一化与偏置:自相关计算结果准不准,全看这两个细节
很多人在用xcorr的时候只关心输出数组,不关心第四个参数scaleopt,结果算出来的自相关数值忽大忽小,或者不同延迟处的方差表现异常。这里把归一化和偏置问题彻底讲透。
3.1 scaleopt 的四档设置:none、biased、unbiased、coeff
xcorr的scaleopt参数一共有四个选项,用一张表说清楚它们都做了什么:
| 选项 | 计算方法 | 适用场景 |
|---|---|---|
'none' | 原始累加和,无缩放 | 主要用于配合功率谱计算,结果包含信号能量信息 |
'biased' | 除以 N | 保证估计一致性,数学性质好,但延迟越大有效样本越少,尾部会出现衰减 |
'unbiased' | 除以 N- | k |
'coeff' | 除以 R(0),使零延迟处为 1 | 最适合观察相关性强弱和周期结构,数值有界,直观 |
对于平稳随机信号,'biased'是理论上的首选,因为它能保证估计值渐进无偏且均方一致。但实际问题中,信号不是无限长的,延迟 k 越接近 N,参与计算的有效样本数就越少,'unbiased'虽然通过除以N-|k|做了补偿,结果却会把尾部那些样本量不足的相关值放大,看起来像出现了很大的波动,其实是噪声被放大了。
3.2 用 'coeff' 还是 'biased':一个实际判断准则
我在做周期检测时,通常先看信号有没有明显直流分量,有则先去均值,然后直接上'coeff'。为什么不用'biased'?因为'coeff'把零延迟归一化到 1,信号在其它延迟处的自相关值都落在 [-1, 1] 之间,我可以直接设一个阈值,比如 0.5,超过阈值的延迟位置基本就能确认是周期峰。这个做法的好处是不用关心信号幅值的绝对大小,只需要看相关性结构。
但如果你需要估计信号的总功率,或者要从自相关函数推导功率谱密度,就不能用'coeff'了,因为它把R(0)归一化了,能量信息被抹掉了。这时候应该用'biased',它的R(0)就是信号的均方值,能直接用于后续谱估计。
3.3 一段验证偏置效应的代码
为了直观展示'unbiased'尾部噪声放大问题,可以做一个快速验证:
N = 1000; x = randn(1, N); % 白噪声 [rb, lags] = xcorr(x, 'biased'); [ru, ~] = xcorr(x, 'unbiased'); subplot(2,1,1); plot(lags, rb); title('biased'); xlim([-50 50]); subplot(2,1,2); plot(lags, ru); title('unbiased'); xlim([-50 50]);白噪声的理想自相关在零延迟处为 1,其余位置为 0。但用'unbiased'处理时,远离零延迟的位置会出现明显的随机起伏,这些起伏并不是信号本身的特征,而是因为延迟越大,用于平均的样本数越少,估计方差越大。这个现象在短信号上尤其明显,所以如果你的数据长度只有几百个点,强烈建议不要用'unbiased'。
4. 三种高频应用场景及完整代码:周期检测、时延估计、噪声鉴别
理论说了一大堆,最终要落到实际场景。自相关和互相关在信号处理里最常见的三个任务就是:从含噪信号中找周期、用互相关计算两个信号的时延、判断信号中的噪声类型。下面逐个给出可复用的代码和结果解读方法。
4.1 周期检测:正弦信号叠加随机噪声,如何用自相关恢复周期
这是自相关最经典的应用。假设你有一段转速传感器信号,里面有一个周期成分和大量噪声,直接看波形可能完全看不出周期,但自相关能把周期成分凸显出来。
fs = 1000; t = (0:4999) / fs; f0 = 20; % 待检测频率 x = sin(2*pi*f0*t) + 1.5*randn(1, 5000); % 去均值后算自相关 x_centered = x - mean(x); [r, lags] = xcorr(x_centered, 'coeff'); % 只取正延迟部分 pos_idx = lags >= 0; r_pos = r(pos_idx); lags_pos = lags(pos_idx); % 寻找除零延迟外的最大峰值 % 先从延迟为0之后找 start_idx = 10; % 跳过前面几个点,避免和零延迟峰混淆 [pks, locs] = findpeaks(r_pos(start_idx:end), lags_pos(start_idx:end)); [~, idx] = max(pks); period_samples = locs(idx); detected_freq = fs / period_samples; fprintf('估计频率: %.2f Hz, 周期: %d 个采样点\n', detected_freq, period_samples);这里要注意,findpeaks找到的第一个主峰位置就是周期对应的延迟。实测时,如果噪声很强,自相关峰值可能出现多个候选位置,处理办法是取前几个显著峰,然后计算它们之间的间隔,再取平均,这样得到的结果更稳。
4.2 时延估计:两个麦克风接收同一信号,用互相关求到达时间差
多麦克风阵列、雷达回波、超声测距,这些任务本质都是时延估计。核心思想是:两个信号是同一源信号的不同延迟版本,互相关函数在真实时延处会出现峰值。
fs = 8000; t = (0:3999) / fs; source = sin(2*pi*300*t) + 0.1*randn(1, 4000); true_delay_samples = 200; x1 = source; x2 = [zeros(1, true_delay_samples), source(1:end-true_delay_samples)]; % 互相关,限制最大延迟范围 maxlag = 500; [r, lags] = xcorr(x2, x1, maxlag, 'coeff'); % 找峰值 [~, idx] = max(abs(r)); estimated_delay = lags(idx); fprintf('真实时延: %d 样本, 估计时延: %d 样本\n', true_delay_samples, estimated_delay);这里要求abs(r)而不是看r本身,是因为如果两个信号反相,互相关峰值会出现在负方向。推荐先确认信号的极性关系,再决定是否加绝对值。另外,maxlag必须设置得比预期的真实延迟大一些,但也不能太大,否则容易把相关峰的定位干扰到别的局部最大值上。
4.3 噪声鉴别:从自相关形状判断信号是白噪声还是有色噪声
白噪声的自相关是一个只在零延迟处尖锐的冲激,而低通滤波后的有色噪声,其自相关在零延迟附近会有一个缓慢衰减的过程。利用这个特征可以快速判断噪声类型,在系统辨识和故障诊断中很实用。
N = 2000; white_noise = randn(1, N); % 构造有色噪声:简单的低通滤波 b = ones(1, 10) / 10; colored_noise = filter(b, 1, white_noise); [rw, lags_w] = xcorr(white_noise - mean(white_noise), 'coeff'); [rc, lags_c] = xcorr(colored_noise - mean(colored_noise), 'coeff'); figure; subplot(2,1,1); plot(lags_w, rw); xlim([-50 50]); title('白噪声自相关'); subplot(2,1,2); plot(lags_c, rc); xlim([-50 50]); title('有色噪声自相关');白噪声那张图,除零延迟外几乎为 0;有色噪声那张图,零延迟两侧有明显的"山包",宽度跟滤波器带宽有关。这个差异肉眼可见,判断起来非常直接。
5. 手写一个自相关函数:从原理到代码,彻底摆脱黑盒
虽然 MATLAB 内置函数很好用,但如果你想加深理解,或者要移植到别的语言,手写一遍自相关是绕不开的。这里给出一个从定义出发的实现,并对比它与xcorr的差异。
5.1 朴素实现的直观写法
根据自相关定义R(k) = Σ x(t)·x(t+k) / (N - |k|),可以用两层循环实现,虽然慢,但能清楚看到每一步在算什么:
function r = my_autocorr(x, maxlag) N = length(x); x = x(:); % 转为列向量 r = zeros(2*maxlag + 1, 1); idx = 1; for k = -maxlag:maxlag % 对齐信号 if k >= 0 x1 = x(1:N-k); x2 = x(1+k:N); else x1 = x(1-k:N); x2 = x(1:N+k); end % 注意这里做了无偏归一化 r(idx) = sum(x1 .* x2) / (N - abs(k)); idx = idx + 1; end end调用方式和内置函数一致,输出长度是2*maxlag + 1。这个实现等价于xcorr(x, 'unbiased'),因为每一次延迟都用实际重叠的样本数做了归一化。它的问题在于复杂度是 O(N·maxlag),信号一长就非常慢,所以只适合教学验证或短序列分析。
5.2 用 FFT 加速的原理与实现
工程上更实用的是用快速傅里叶变换把时域卷积变成频域相乘。自相关与功率谱的对应关系——维纳-辛钦定理——告诉我们:自相关函数的傅里叶变换等于功率谱密度。反过来,先计算信号的功率谱,再逆傅里叶变换,就能得到自相关函数。这个做法的复杂度是 O(N log N),对长信号几乎是唯一可行的方案。
N = 10000; x = randn(1, N); x = x - mean(x); % 为使线性相关而非循环相关,需要补零到 2N-1 nfft = 2^nextpow2(2*N - 1); X = fft(x, nfft); Sxx = X .* conj(X); % 功率谱 r_ifft = ifft(Sxx); % 截取实际相关部分,并做无偏归一化 r_ifft = r_ifft(1:N); lags = 0:N-1; r_unbiased = r_ifft ./ (N - lags); % 比较与内置 xcorr 的结果 [r_builtin, lags_builtin] = xcorr(x, 'unbiased'); builtin_positive = r_builtin(lags_builtin >= 0); max_diff = max(abs(builtin_positive - r_unbiased)); fprintf('与内置函数的最大偏差: %.2e\n', max_diff);注意这里两个关键点:一是补零到2N-1以上,否则得到的是循环相关而不是线性相关;二是除法处理无偏归一化时,分母是N - lag,也就是每个延迟对应的有效样本数。实测最大偏差通常在1e-12量级,完全来自浮点运算误差。
5.3 手写时最容易踩的两个坑
第一个坑是忘记去均值。如果你在定义中使用中心化自相关(即自协方差),却不对信号做去均值处理,结果会整体上移,零延迟处尤其明显。第二个坑是补零长度不足。FFT 方法要求 nfft 至少是2N-1,否则时间混叠会把尾部错误地折叠到头部,导致相关函数左右不对称。这两个坑我在早期移植代码时都踩过,定位问题花了不少时间,这里提前帮你排掉。
6. 从 xcorr 结果到功率谱密度:自相关的一个重要应用
相关分析并不止步于算一条曲线。前面反复提到维纳-辛钦定理,这里实际做一遍,把自相关和功率谱的换算链路打通。这在信号处理中非常常用:估算信号的功率谱,可以直接 FFT 后取模方,也可以通过自相关的傅里叶变换来估计,后者在一些现代谱估计方法中效果更好,比如 Welch 法和周期图法其实都是对自相关思路的改进。
6.1 用自相关法估计功率谱的示例
fs = 1000; t = (0:4095) / fs; x = sin(2*pi*100*t) + 0.5*sin(2*pi*250*t) + randn(1, 4096); % 计算有偏自相关,保证非负定性 [r, lags] = xcorr(x, 'biased'); % 对自相关做 FFT 得功率谱 Pxx = fft(fftshift(r)); Pxx = abs(Pxx(1:length(Pxx)/2+1)); freq = linspace(0, fs/2, length(Pxx)); % 参考:直接周期图 Pxx_direct = pwelch(x, 256, 128, 256, fs); plot(freq, 10*log10(Pxx/Pxx(1)), 'linewidth', 1.5); hold on; plot(linspace(0, fs/2, length(Pxx_direct)), 10*log10(Pxx_direct), 'r');这里的小技巧是,xcorr的输出默认峰值在数组中间,要先用fftshift把零延迟移到数组开头,再做 FFT,得到的频率轴才正确。自相关法得到的谱曲线比周期图平滑,但分辨率稍低;周期图则恰好相反。二者互补,实际工程中常联合使用。
6.2 为什么自相关法在某些场景下更稳
直接周期图法有个问题:对信号加窗后,频谱泄漏和窗函数旁瓣会对弱信号成分造成干扰,而且数据越长,谱估计的方差并不收敛,还是一样大。而把自相关作为中间步骤,相当于先对信号做了统计平均,再进傅里叶变换,谱估计的方差相对可控。这也是现代谱估计中的基本直觉。
不过自相关法也不是万能的。当信号包含强周期分量时,自相关的旁瓣会产生频谱泄漏,导致弱信号被埋没。这时候可以考虑先做频域平滑,或者在自相关上加窗(比如取前 M 个延迟再做 FFT),这些都是后话,但值得知道自相关这条路不是越走越长就好,而是越"精炼"越稳。
7. 实操中的若干坑:关于数据长度、直流分量和双变量扩展
这部分单纯是经验总结,没有任何理论推导,全是教训。我在不同项目里反复和自相关打交道,遇到过的问题可以归为这几类。
7.1 数据长度不足导致周期性误判
自相关对数据长度极其敏感。假设信号周期是 100 个采样点,你只有 150 个点的数据,那么自相关函数里能观察到的峰非常有限,而且第 2 个峰只有 50 个样本参与计算,幅度和形状都会变形,容易被误判为噪声。经验法则是:数据长度至少要有待检测周期的 5 到 10 倍,否则宁可先用带通滤波把信号处理干净,再做相关分析。
7.2 直流分量对互相关时延估计的干扰
进行两个信号的互相关时,如果两路信号都带有不同的直流分量,峰值位置虽然一般不受影响,但相关曲线上会叠加一个缓慢变化的斜坡,严重时会让峰值定位变得模糊。实测中,先分别对两路信号做去均值再算互相关,是最稳妥的做法。我在做超声回波时延估计时,第一版代码忘了去均值,结果峰值偏移了 3 个采样点,换算成距离误差大约 0.5 毫米,后来加上去均值后就完全正常了。
7.3 从单变量自相关到双变量空间自相关:思路的扩展
热搜词里出现了"双变量空间自相关",这个方向也值得提一句。地理信息或空间统计中的双变量自相关,本质上是把时间序列的自相关推广到空间域和时间-空间联合域。比如分析两个变量的空间分布是否相关,可以用双变量 Moran's I,它在计算矩阵形式上和互相关有相通之处,但多了空间权重矩阵的概念。如果你已经熟悉了 MATLAB 的互相关计算,理解空间自相关会有天然的优势,因为核心思路都是"衡量一组变量在不同偏移下的相关结构",只是偏移的维度从时间轴换成了空间邻接关系。在 MATLAB 中处理空间自相关,可以用corrcoef配合空间权重矩阵自行实现,或者通过 Statistics and Machine Learning Toolbox 完成,这里不展开,但值得作为一个进阶方向思考。
7.4 长序列计算的效率问题
当信号长度达到几十万点以上,直接调xcorr虽然有内置优化,但内存开销也不小。实测一段 100 万点的信号,xcorr输出的数组大约有 200 万个元素,占用约 16 MB 内存,这还不算中间变量。如果还要做不同maxlag的多次尝试,建议尽量限制maxlag,或者用 5.2 节中的 FFT 思路自行实现。特别是做实时处理时,控制maxlag往往比优化代码本身更有效。
8. 收尾:关于"相关函数"这个压缩包的几点忠告
如果你手头也刚好拿到这样一个命名杂乱、素材零散的压缩包,我的建议是不要急着运行任何.m文件,先自查三个问题:第一,压缩包里是否有说明文件或注释,标明了计算用的是哪种归一化方式;第二,涉及的信号是随机信号还是确定信号,是否满足平稳性假设;第三,计算自相关之前,有没有做预处理,比如去均值、去趋势、滤除直流。
这三个问题在相关函数分析里是决定结果可靠性的关键,比代码本身重要得多。我遇到过太多次"明明代码一模一样,结果却对不上"的求助帖,最后排查出的原因往往不是算法问题,而是前置条件没有对齐。
回到 MATLAB 本身,自相关和协方差的计算已经足够成熟,xcorr、xcov、autocorr这些函数在各种场景下覆盖了绝大多数需求。我的建议是:概念上用统计学定义去理解,工程上用xcorr的'coeff'参数去落地,需要谱分析时切换到'biased',然后结合findpeaks或pwelch做后续处理。这套组合拳基本能覆盖周期检测、时延估计、噪声鉴别、谱估计等常见任务,也是我这些年使用频率最高的一套流程。希望对你有帮助。
本文还有配套的精品资源,点击获取