最近在相关高斯噪声下的时延估计，，进行仿真实验，需要生成两路相关的高斯白，由于从来没接触过，在MATLAB官网找到一个函数，如下。。求教，查了很多论文，发现直接互相关法在相关噪声背景下几乎是失效的，可是我按照这个函数写的时延估计程序，互相关的方法还能得出很好的结果。另外，做基于四阶累积量与RLS自适应滤波结合估计时延时，误差系数曲线异常，找了很久的资料，也还没解决，，，求各位前辈指教！！！
（Multisynchrosqueezing Transform - File Exchange - MATLAB Central https://ww2.mathworks.cn/matlabcentral/fileexchange/68571-multisynchrosqueezing-transform）
function [result, sampRpp] = correlatedGaussianNoise(Rpp, nSamp)
%% generates correlated 0-mean Gaussian vector process
%生成相关的0均值高斯向量过程
% inputs: Rpp - correlation matrix. must be positive definite %相关矩阵。必须是正定的 大小决定输出向量
% and size determines output vector
% nSamp - number of independent samples of correlated process相关过程的独立样本数量
%
% output: result - matrix of dimension Rpp rows X nSamp cols. 矩阵的尺寸Rpp行X nSamp cols
%
% result has form:
% < ------ nSamp ------->
% ------------------------ data is correlated along 数据是相关的
% | | | all rows for each column
% | | output |
% p | data matrix | data is independent along
% | | | all columns for a given row给定行的所有列
% | | |
% < -------- nSamp ------>
%
% example: use following correlation matrix (output implicitly 3 rows) 示例:使用以下相关矩阵(隐式输出3行)
% [1 0.2 0.2]
% [0.2 1 0.2]
% [0.2 0.2 1 ]
%
% and 1e5 samples to check function
%
%%
% n = 3; Rpp = repmat(0.2, [n,n]); Rpp(1:(n+1):n^2) = 1;
% disp(Rpp)
% nSamp = 1e5;
% [x, sampR] = correlatedGaussianNoise(Rpp, nSamp);
% disp(sampR)
%
%% -----------------------------------------------------
% michaelB
%
%% algorithm
% check dimenions - add other checking as necessary...
if(ndims(Rpp) ~= 2),
result = [];
error(\'Rpp must be a real-valued symmetric matrix\');
return;
end
% symmeterize the correlation matrix
Rpp = 0.5 .*(double(Rpp) + double(Rpp\'));
% eigen decomposition
[V,D] = eig(Rpp);
% check for positive definiteness
if(any(diag(D) <= 0)),
result = [];
error(\'Rpp must be a positive definite\');
return;
end
% form correlating filter
W = V*sqrt(D);
% form white noise dataset
n = randn(size(Rpp,1), nSamp);
% correlated (colored) noise
result = W * n;
% calculate the sample correlation matrix if necessary
if(nargout == 2),
sampRpp = 0;
for k1=1:nSamp,
sampRpp = sampRpp + result(:,k1) * (result(:,k1)\');
end
% unbiased estimate
sampRpp = sampRpp ./ (nSamp - 1);
end
% % % % 基于四阶累积量与RLS自适应滤波的时延估计方法（FOC-RLS）%%%%%% clc clear all close all N = 500; fs=1000; a=-0.5; nn = 2; Rpp = repmat(0.9, [nn,nn]); Rpp(1:(nn+1):nn^2) = 1; disp(Rpp); nSamp=500; [result, sampRpp] = correlatedGaussianNoise(Rpp, nSamp); n1=0.4*result(1,:); %result为2行矩阵，第一行 为噪声1，第二行 为噪声2 n2=0.4*result(2,:);%%系数调整信噪比的大小 t=linspace(1,10,N); sin1=exp(a*t); sin2=zeros(1,length(sin1)); delay= 100; sin2((delay+1):length(sin1))=sin1(1:length(sin1)-delay); nn1=sin(2*pi*10*t); nn2=zeros(1,length(nn1)); nn2((delay+1):length(nn1))=nn1(1:length(nn1)-delay); single_1=sin1+n1+nn1; %single_A single_2=sin2+n2+nn2;%%single_B snr1=SNR_singlech(sin1,single_1); snr2=SNR_singlech(sin2,single_2); figure, plot(t,single_1,\'color\',\'b\'); hold on; plot(t,single_2,\'color\',\'g\'); legend(\' 信号A\',\'信号B\'); X1=fft(single_1,2*N-1); X1J=conj(X1); X2=fft(single_2,2*N-1); X2J=conj(X2); X12=X2.*X1J; R12=fftshift(real(ifft(X12))); r12max=max(abs(R12)); R12=R12 / r12max; lags=-499:499; figure, plot(lags,R12,\'color\',\'b\');title(\'基本互相关\'); xCum4 = cum4est (single_1, 90, 20, 0, \'biased\', 10, 10);%自四阶累积 xCum4Size = size(xCum4); %返回数组维度 absmax1=max(abs(xCum4)); xCum4=xCum4/absmax1; xCum4_hu = cum4est (ifft(X12), 90, 20, 0, \'biased\', 10,10);%互四阶累积 xCum4Size_hu = size(xCum4_hu); %返回数组维度 absmax=max(abs(xCum4_hu)); xCum4_hu=xCum4_hu/absmax; figure, subplot(2,1,1) plot(xCum4,\'color\',\'b\') ; title(\'信号1四阶自累积量运算\'); xlabel(\'子频带数\') subplot(2,1,2) plot(xCum4_hu,\'color\',\'b\') ; title(\'互四阶累积量运算\'); xlabel(\'子频带数\') jieguo=xcorr(xCum4,xCum4_hu,\'unbiased\'); %%频域的累积量互相关 fore=ifft(jieguo); %%%反变换到时域 figure, plot(jieguo,\'color\',\'b\') ; title(\'四阶累积量互相关\'); [v1,location1]=max(abs(fore)); delaytime_foc=(fix(length(jieguo)/2)-location1)/fs %经过累积计算换算出信号的延迟时间 [v2,location2]=max(abs(R12)); location2=lags(location2); delaytime_cc=location2/fs %互相关换算出信号的延迟时间 % % % % % % % % RLS % % % % % % % % % % % % % % % % x1=xCum4; %产生高斯随机系列 d1=sin1; % 期望输出 [y1 e1 w1] = RLS(x1, d1, 32); figure, subplot(2,1,1) plot(e1,\'color\',\'b\') ; title(\'y1_error\'); subplot(2,1,2) plot(w1(end,:),\'color\',\'b\') ; title(\'y1 滤波器权系数\'); x2=y1\'-xCum4_hu; d2=sin2; [y2 e2 w2] = RLS(x2, d2, 32); figure, subplot(2,1,1) plot(e2,\'color\',\'b\') ; title(\'误差曲线\');xlabel(\'迭代次数\') ; ylabel(\'输出误差估计\') ; subplot(2,1,2) plot(w2(end,:),\'color\',\'b\') ; title(\'滤波器权系数\'); xlabel(\'迭代次数\') ; ylabel(\'权矢量的值\') ; [v3,location3]=max(abs(e2)); delaytime_RLS=(location3)/fs %经过累积计算换算出信号的延迟时间