相關噪聲下時延估計

最近在相關高斯噪聲下的時延估計，，進行仿真實驗，須要生成兩路相關的高斯白，因爲歷來沒接觸過，在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 %通過累積計算換算出信號的延遲時間