% Autocorrelation function detection
clear all
close all
fs=10000;
N=1000;
x=0:1/fs:(fs-1)/fs;
y0=sin(2*pi*1000*x);
n=fs;
x1=normrnd(0,1,1,n);
%散粒噪声
lamdac=2;Tmaxc=1000;
ic=1;Tc(1)=random('exponential',lamdac);
while(Tc(ic)<Tmaxc)
Tc(ic+1)=Tc(ic)+random('exponential',lamdac);
ic=ic+1;
end
Tc(ic)=Tmaxc;xc=0:1:ic;wc(1)=0;
for pc=1:ic
wc(pc+1)=Tc(pc);
end
tc=(1/fs:1/fs:1)*Tmaxc;
ac=0.1;
Wc=fs/Tmaxc*wc;
Wc=floor(Wc);
yc=zeros(length(Wc)-1,fs);
for k1c=2:1:length(Wc);
kc=Wc(k1c);
t1c=1:1:fs;
ypc=exp(-(t1c-kc)*ac);
for kpc=1:1:kc-1
ypc(kpc)=0;
end
yc(k1c-1,:)=ypc;
end
x2=zeros(1,fs);
for k2c=1:1:fs
x2(1,k2c)=sum(yc(:,k2c));
end
% 高斯色噪声
p=0.5;
f0=0.05;
a1=-2*p*cos(2*pi*f0);
a2=p^2;
x3=zeros(1,n);
for i=3:1:n;
x3(i)=a1*x3(i-1)-a2*x3(i-2)+x1(i);
end
% 混合噪声形成
y=y0+x1;
yfft=abs(fftshift(fft(y)));
f=(-length(yfft)/2:(length(yfft)/2-1))*fs/length(yfft);
figure(1)
plot(x,y)
ylabel('幅值');
xlabel('时间');
title('混合信号时域波形');
figure(2)
plot(f,yfft)
ylabel('yfft');
xlabel('f');
title('混合信号频谱');
% 计算输入信噪比
snr1=0;
ps1=sum(sum((y0-mean(mean(y0))).^2));
pn1=sum(sum((y-y0).^2));
snr1=abs(10*log10(ps1/pn1));
ys=y;
% 一重自相关
[y1,lags]=xcorr(ys,'unbiased');
t=((-length(lags))/2/length(lags):1/length(lags):(length(lags)/2-1)/length(lags))*2;
ys1=2*y1;
figure(6)
plot(t,ys1)
ylabel('幅值');
xlabel('时间');
title('自相关输出');
% 计算自相关输出信噪比
fs1=2*fs;
ts1=0:1/fs1:(fs1-2)/fs1;
yr1=sin(2*pi*2*1000*ts1);
snr2=0;
ps2=sum(sum((yr1(5000:15000)-mean(mean(yr1(5000:15000)))).^2));
pn2=sum(sum((ys1(5000:15000)-yr1(5000:15000).^2)));
snr2=abs(10*log10(ps2/pn2));
% 二重自相关
[y2,lags]=xcorr(ys1,'unbiased');
t1=((-length(lags))/2/length(lags):1/length(lags):(length(lags)/2-1)/length(lags))*4;
ys2=2.2*y2;
figure(7)
plot(t1,ys2)
ylabel('幅值');
xlabel('时间');
title('二重自相关输出');
% 计算二重相关输出信噪比
fs2=4*fs;
ts2=0:1/fs2:(fs2-4)/fs2;
yr2=sin(2*pi*1000*ts2);
snr3=0;
ps3=sum(sum((yr2(10000:20000)-mean(mean(yr2(10000:20000)))).^2));
pn3=sum(sum((ys2(10000:20000)-yr2(10000:20000))).^2);
snr3=abs(10*log10(ps3/pn3));
