相关函数

% 功能：
%       1.导入txt格式的两信号，
%       2.求其各自的自相关函数Rxx，自协方差函数Cxx，自相关系数/归一化自协方差函数
%       3.求两信号的互相关函数Rxy，互协方差函数Cxy,，互相关系数/归一化互协方差函数，延时tao
%       4.求信号的功率（求均方值，或者自相关函数(延时tao=0),自功率谱的IFFT(延时tao=0，注意区别双边谱和单边谱)）
% 各物理意义如下：
% 自相关函数Rxx：反映了噪声x(t)在不同时刻t1,t2取值的相关程度
%               延时为0时的R值（最大值）= 均方值，反映了随机噪声的平均功率（直流分量功率+交流分量功率）
%               延时为无穷大时的R值 = 均值^2, 反映了随机噪声的直流分量的功率
% 归一化相关函数：也反映了噪声x(t)在不同时刻t1,t2取值的相关程度，
%                但消除了随机噪声的幅度和功率的影响，更准确反映相关程度
% 归一化自相关函数 = Rxx/Rxx(0) = Rxx/max(Rxx) = Rxx/ ((1/length(x))*sum(x.^2)); 
%                  ='coeff'
% 归一化互相关函数 =  Rxy/sqrt(Rxx(0).*Ryy(0)) = Rxy/sqrt(max(Rxx).*max(Ryy(0)))
%                  ='coeff'
% 注意：均方值因为除了个N，所以对应的是biased，注意对应关系
% 自协方差函数Cxx ：反映了同一随机噪声x(t)的交流分量在不同时刻t1,t2取值的相关程度
%                Cxx  = Rxx-(mean(x))^2
% 互协方差函数Cxy：反映了两随机噪声x(t)的交流分量在不同时刻t1,t2取值的相关程度
%                Cxx  = Rxy-mean(x)*mean(y)
% 归一化自协方差函数Cxx_nom = Cxx./var(x)
% 归一化互协方差函数Cxy_nom = Cxy./(sqrt(var(x)).*sqrt(var(y)))
% 编辑者：lily
% 日期：2019,4,13
clc;
clear;
close all; 
% =========== input signal ==========================
Fs = 12000;
x = load('err mic1 - off.txt');
y = load('err mic2 - off.txt');
% =============== AutoCorr/CrossCorr ===============================
% N ：累加平均次数
% m ：延时序号
% 'biased'：有偏估计，将'none'/(N)
% 'unbiased'：无偏估计，将'none'/(N-|m|)
% 'coeff'：将'none' 作归一化处理,除以最大值。
%          Rxx_max: m = 0时的值，也是均方值                 (  Rxx_max = Rxx(0) )
%          Rxy_max: 两信号m = 0的自相关函数相乘，再开根号。（  Rxy_max = sqrt(Rxx(0)*Ryy(0)) ）
% 'none'：不作归一化处理
%  biased虽是有偏估计，但渐近一致估计，估计量的方差小于无偏估计，实际中多用有偏估计
%  归一化相关函数 = 用'none'先计算，再求归一化 = 'coeff'
str = {
    'biased','unbiased','coeff','none'};
ind = 1;
[Rxy, c] = xcorr(x, y, str{
    ind});
[Ryy, d] = xcorr(y, y, str{
    ind});
[Rxx, b] = xcorr(x, x, str{
    ind});
tao = b/Fs;       %tao: 偏移的时间，用偏移个数乘采样间隔
% =============== figure ===============================
figure;
plot(tao, Rxx);
title(str{
    ind});
grid;
figure;
plot(tao, Rxy);
title(str{
    ind})
grid;
% ===============求延时 =====================================
 [~,indexMax] = max(Rxx);
 disp(['Time lag = ' num2str(tao(indexMax))])
 % ============= 归一化自相关函数rho_xx ==============================
rho_xx = Rxx./max(Rxx) ;
rho_xxx = Rxx./ ((1/length(x))*sum(x.^2));
figure;
plot(rho_xx);
hold on;
plot(rho_xxx);
legend('用最大值算','用均方值算');
title('归一化自相关函数rho_xx ');
% =============== 归一化互相关函数rho_xy =================================
rho_xy = Rxy/sqrt(max(Rxx).*max(Ryy));
figure;
plot(rho_xy);
title('归一化互相关函数rho_xy ');
 % =============== 自协方差函数Cxx ===============================
Cxx = Rxx-(mean(x))^2;
figure;
plot(tao, Cxx);
title('自协方差函数Cxx')
grid;
 % =============== 互协方差函数Cxy ===============================
Cxy = Rxy-mean(x)*mean(y);
figure;
plot(tao, Cxy);
title('互协方差函数Cxy ')
grid;
% =========== 归一化自协方差函数Cxx_nom ============================
Cxx_nom = Cxx./var(x);
figure;
plot(tao,Cxx_nom);
title('归一化自协方差函数Cxx_nom')
grid;
% ============ 归一化互协方差函数Cxy_nom=============================
Cxy_nom = Cxy./(sqrt(var(x)).*sqrt(var(y)));
figure;
plot(tao,Cxy_nom);
title('归一化互协方差函数Cxy_nom')
grid;

已知平稳随机过程的功率谱，求自相关函数

%  功能：已知平稳随机过程的功率谱，求自相关函数
%  平稳随机过程的功率谱连续时间表达形式：Fw = 2*pi* dirac(w);
%  验证频域为1的ifourier是单位冲击函数。
% 编辑者：lily
% 日期：2019,4,18
clear;
clc;
close all;
% ======================= input signal ==========================
syms w t  %定义一个符号变量x, 后续可以做一些符号操作
Fw = 2*pi* dirac(w);
ft = ifourier(Fw,w,t);
pretty(ft);  % 表达式的形式