%%%********对于不同的加权函数下 不同的信噪比对时延的影响*********
clear all; clc; close all;
N=1024; %信号长度
Fs=500; %采样频率
n=0:N-1;
t=n/Fs; %时间序列
a1=2; %信号幅度
a2=2;
d=3; %延迟点数
x1=a1*cos(2*pi*10*n/Fs); %信号1
%x1=x1+randn(size(x1)); %加随机噪声--[randn(size(x1)) 返回一个和A有同样维数大小的随机数组]
x1=awgn(x1,20); %加信噪比为20的噪声
x2=a2*cos(2*pi*10*(n+d)/Fs); %信号2
%x2=x2+randn(size(x2)); %加随机噪声
x2=awgn(x2,20); %加信噪比为20的噪声
subplot(211); %设置子图
plot(t,x1,'r'); %作图 红色线条
hold on; %在新画图像之后不想覆盖原图像就要加上hold on这句话
plot(t,x2,':');
%legend('x1信号', 'x2信号'); % 设置图例的线条和patches的函数，线条和字幕的设置等
xlabel('t/s');ylabel('x1(t) /x2(t)');
title('初始信号');
grid on; %添加网格线
hold off
%互相关函数
X1=fft(x1,2*N-1); %傅里叶变换
X2=fft(x2,2*N-1); %FFT
%Cxy=fftshift(ifft(Sxy)); %求互相关函数ifft(逆)，fftshift逆变换后矩阵左右转换
%Cxy=fftshift(real(ifft(Sxy)));% real取复数实部
%subplot(212);
%f=n*Fs/N-Fs/2;f2=n*Fs/N;mag=abs(Cxy);stem(f,mag);
%t1=(0:2*N-2)/Fs; %从0开始的，x1是2047个数，所以是 2*N-2 时间取数的原理
%%%%*************** PHAT-GCC 加权*********************************
%Sxy=X1.*conj(X2); %共轭相乘得互相关函数Sxy
%Cxy=fftshift(ifft(Sxy./abs(Sxy)));
%subplot(212);
%t1=(-N+1:N-1)/Fs;
%5%%*************************************************************
%%%*****************ROTH-GCC 加权*********************************
%X12=X1.*conj(X2);
%X11=X1.*conj(X1);
%Cxy=real(fftshift(ifft(X12./abs(X11))));
%t1=(-N+1:N-1)/Fs;
%%%*************************************************************
%%%*****************SCOT-GCC 加权*********************************
X=X1.*conj(X2);
X11=X1.*conj(X1);
X22=X2.*conj(X2);
Y=sqrt(X11.*X22);
Cxy=real(fftshift(ifft(X./Y)));
t1=(-N+1:N-1)/Fs;
%%%*************************************************************
plot(t1,Cxy,'g');
xlabel('时间/t');ylabel('Cxy(t)');grid on
%gtext('SNR=-20'); 波线加标注
[max,location]=max(Cxy);%求出最大值max,及最大值所在的位置（第几行）location;
if location>0
dd=location-1;
end
if location>(2*N-1)/2
dd=location-(2*N-1)/2;
end
Delay=dd/Fs
%%%%************如下算法是相同信噪比不同算法的比较***********
clear all; clc; close all;
b0=[];
c0=[];
for snr=-21:3:21 %不同信噪比计算不同的方差，循环
a0=[];
for i=1:50 %循环程序多少遍计算出的一个方差
N=1024; %长度
Fs=500; %采样频率
n=0:N-1;
t=n/Fs; %时间序列
a1=2; %信号幅度
a2=2;
d=3;
x1=a1*cos(2*pi*10*n/Fs); %信号1
x1=awgn(x1,snr); %加噪声
%x1=x1.*hamming(max(size(x1)))';%加窗
x2=a2*cos(2*pi*10*(n+d)/Fs); %信号2
x2=awgn(x2,snr);
%x2=x2.*hamming(max(size(x2)))';%加窗
subplot(511);
plot(t,x1,'r');
hold on;
plot(t,x2,':');
legend('x1信号', 'x2信号');
xlabel('时间/s');ylabel('x1(t) x2(t)');
title('原始信号');grid on;
%互相关函数
tic
X1=fft(x1,2*N-1);
X2=fft(x2,2*N-1);
Sxy=X1.*conj(X2);
%Cxy=fftshift(ifft(Sxy));
%%%%% GCC ****************
Cxy=fftshift(ifft(Sxy));
subplot(512);
t1=(-N+1:N-1)/Fs;
plot(t1,Cxy,'b');
title('GCC');xlabel('t/s');ylabel('Cxy');grid on;
[max1,location1]=max(Cxy);
d1=location1-N
Delay1=d1/Fs
%%%%% ************************PHAT-gcc*************************
%X1=fft(x1,2*N-1);
%X2=fft(x2,2*N-1);
%Sxy=X1.*conj(X2);
%Pxy=fftshift(ifft(Sxy./abs(Sxy)));
%subplot(513);
%t1=(-N+1:N-1)/Fs;
%**************************************************************
%%%*****************SCOT-GCC 加权*********************************
X=X1.*conj(X2);
X11=X1.*conj(X1);
X22=X2.*conj(X2);
Y=sqrt(X11.*X22);
Cxy=real(fftshift(ifft(X./Y)));
t1=(-N+1:N-1)/Fs;
%%%*************************************************************
plot(t1,Cxy,'b');
title('PHAT-GCC');xlabel('t/s');ylabel('Pxy');grid on;
[var1,location2]=max(Cxy);
d2=location2-N
Delay2=d2/Fs
A=[a0,Delay2];
a0=A; % 每个delay输出一个一维数组
end
VX=var(A) %在不同的信噪比下得出不同的方差，求出信噪比和方差的关系图
%fprintf('方差是：%d\n',VX)
B=[b0,VX];
b0=B;
C=[c0,snr];
c0=C;
%fprintf('SNR是：%d\n',snr)
end
figure %%%根据计算出的两组数据，SNR与信噪比关系作图
plot(C,B,':kh')
xlabel('SNR(db)');ylabel('方差');title('SCOT-GCC')
set(gca,'Xtick',[-21:3:21]);
set(gca,'Ytick',[0.01:0.1:2]);
