MATLAB之相关函数

谁践踏了优雅 2022-02-20 14:06 701阅读 0赞

### 文章目录 ###

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

# 相关函数 #

% 功能：
    %       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);  % 表达式的形式