clear all
clc
format long
f0=500;
c=1500;
d=0.001*c/f0;
Fs=10000;
lth=1024;               % length of FFT
sita=30;
SNR=-70;
count=16;    
t0=d/c*sin(sita*pi/180);
S=zeros(count,lth);  

% ********產生信號********

for k=1:3
    for j=1:lth
        s(k,j)=exp(i*(2*pi*f0*(j-1)/Fs-2*(k-1)*pi*f0*t0));
    end  
end
subplot(4,2,1)
j=1:lth;
t=(j-1)/Fs;
plot(t,s(1,:));
axis([0 0.025 -1.5 1.5])

 
 %*******根據信噪比產生噪聲********
    
 for k=1:count
      n(k,:)=randn(1,length(s(1,:)));         
 end
 for k=1:count;
     A(k,:)=sqrt(10^(-SNR/10)*sum(abs(s(1,:)).^2)/sum(abs(n(k,:)).^2));
 end
 n1=0;
 for k=1:count-2;
     n1= n1+A(k,:)*n(k,:);
end
n2=0;
for k=2:count-1;
     n2=n2+A(k,:)*n(k,:);
end
n3=0;
for k=3:count;
     n3=n3+A(k,:)*n(k,:);
end

p11=s(1,:)+A(1,:)*n(1,:);
p1=(count-2)*s(1,:)+n1;
P1=fft(p1);
N1=fft(n1);
subplot(4,2,3 )
plot(t,p11)
axis([0 0.025 -10000 10000]);
subplot(4,2,2)
[S,f]=pwelch(s(1,:),[],[],1024,Fs);
plot(f,S)
subplot(4,2,4)
[Pxx1,f]=pwelch(p11,[],[],1024,Fs);
plot(f,Pxx1)

p2=(count-2)*s(2,:)+n2;
P2=fft(p2);
N2=fft(n2);
p3=(count-2)*s(3,:)+n3;
P3=fft(p3);
N3=fft(n3);
subplot(4,2,5)
plot(t,p1),axis([0 0.025 -20000 20000]);
subplot(4,2,6)
[Sxx,f]=pwelch(p1,[],[],1024,Fs);
plot(f,Sxx)

w=(1:length(P1))*2*pi*Fs/length(P1);
A=real(P3).*sin(2*w*t0)+imag(P3).*cos(2*w*t0)-imag(P1);
B=real(P2).*sin(w*t0)+imag(P2).*cos(w*t0)-imag(P1);
C=real(P2).*cos(w*t0)-imag(P2).*sin(w*t0)-real(P1);
D=real(P3).*cos(2*w*t0)-imag(P3).*sin(2*w*t0)-real(P1);
Res=((abs(P2).^2-abs(P1).^2-abs(N2).^2+abs(N1).^2).*A-(abs(P3).^2-abs(P1).^2-abs(N3).^2+abs(N1).^2).*B)./(2*(count-2)*(C.*A-D.*B));
Ims=((abs(P3).^2-abs(P1).^2-abs(N3).^2+abs(N1).^2).*C-(abs(P2).^2-abs(P1).^2-abs(N2).^2+abs(N1).^2).*D)./(2*(count-2)*(C.*A-D.*B));
sout=ifft(Res+i*Ims);
subplot(4,2,7)
j=1:lth;
plot((j-1)/Fs,real(sout))
axis([0 0.025 -1.5 1.5])
subplot(4,2,8)
[Sxx1,f]=pwelch(real(sout), [], [],1024,Fs);
plot(f,Sxx1)