%Multiplepath,信道多徑數(shù)
%N，反射點的個數(shù)（通常取為32）
%N0,低頻振蕩器個數(shù):N0=N/4
%fc,載波頻率
%fm,多普勒頻移 
%fs,采樣頻率
%velocity,移動臺移動速度
%framelength,一幀長度
%inphase,瑞利信道的同相分量
%quadrature,瑞利信道的正交分量
%H,瑞利信道系統(tǒng)函數(shù)

%------------------初始化-----------------
Multiplepath = 8;
N = 32; N0=N/4;
fc = 9e8;
velocity = 60; % unit is km/h
fm = (velocity*1000/3600) * fc / 3e8;
fs = 9.6e3;
framelength = 100000;
inphase = zeros(Multiplepath, framelength);
quadrature = inphase;
H = inphase;



%---------------------------------------
for row=0:Multiplepath-1
%fia = rand(1, N0)*2*pi;
%fia1 = rand(1, N0)*2*pi;
fia = random('unif',-pi,pi,1,N0);
fia1 = random('unif',-pi,pi,1,N0);
%pusia = random('unif',-pi,pi,1,N0);
alpha = zeros(1, N0);
for n=0:(N0-1)
alpha(1, n+1) = (2*pi*n-pi+random('unif',-pi,pi))/N;
end

for column=0:framelength-1
tempinphase = zeros(1, N0);
tempquadrature = zeros(1, N0);

for n=1:N0
tempinphase(1, n) = cos(2*pi*fm*cos(alpha(1,n))*(column/fs)+fia(1,n));
tempquadrature(1, n) = cos(2*pi*fm*sin(alpha(1,n))*(column/fs)+fia1(1,n));
end
inphase(row+1, column+1) = sum(tempinphase);
quadrature(row+1, column+1) = sum(tempquadrature);
end
end
H = sqrt(1/N0) * (inphase + sqrt(-1)*quadrature);

%-----------檢驗所得信道的自相關、互相關函數(shù)
%J = zeros(1, 100000);
%for i=0:9999
% J(i+1) = besselj(0, 2*pi*fm*i/fs);
%end

save = sqrt(2/N0)*inphase(1,:);
save1 = sqrt(2/N0)*inphase(1,:);
%test = zeros(1, 1000);
test = zeros(1,20000);
result = test;
for i=1:20000
%for i=1:1000
if i==1
%result(i) = mean(save.*save1);
result(i) = mean(save(1:20000).*save1(1:20000));
%result(i) = sum(save)*sum(save1);
end
if i>1
test = save1(i:i+20000-1);
test = save(1:20000) .* test;
result(i) = mean(test);
%test = [save1(i:1000) save1(1:i-1)];
%test = save.*test;
%result(i) = mean(test);
%test = sum(save)*sum(test);
%result(i) = test;
end
end
%result
%v = [0, 1000, -0.6, 0.6];
v = [0,50000,-0.6,0.6];
plot(1/2*result, 'r'), axis(v)
hold on;
%--------------------下面為比較------------------
for row=0:Multiplepath-1
%fia = rand(1, N0)*2*pi;
%fia1 = rand(1, N0)*2*pi;
fia = random('unif',-pi,pi,1,N0);
fia1 = random('unif',-pi,pi,1,N0);
alpha = zeros(1, N0);
for n=0:(N0-1)
alpha(1, n+1) = (2*pi*n-pi+random('unif',-pi,pi))/N + (2*pi*row-pi+random('unif',-pi,pi))/N;
end

for column=0:framelength-1
tempinphase = zeros(1, N0);
tempquadrature = zeros(1, N0);

for n=1:N0
tempinphase(1, n) = cos(2*pi*fm*cos(alpha(1,n))*(column/fs)+fia(1,n));
tempquadrature(1, n) = sin(2*pi*fm*sin(alpha(1,n))*(column/fs)+fia1(1,n));
end
inphase(row+1, column+1) = sum(tempinphase);
quadrature(row+1, column+1) = sum(tempquadrature);
end
end
H = sqrt(1/N0) * (inphase + sqrt(-1)*quadrature);

%-----------檢驗所得信道的自相關、互相關函數(shù)
%J = zeros(1, 100000);
%for i=0:9999
% J(i+1) = besselj(0, 2*pi*fm*i/fs);
%end

save = sqrt(1/N0)*inphase(4,:);
save1 = sqrt(1/N0)*quadrature(4,:);
test = zeros(1, 50000);
result1 = test;
for i=1:50000
if i==1
result1(i) = mean(save(1:50000).*save1(1:50000));
end
if i>1
test = save1(i:i+50000-1);
test = save(1:50000) .* test;
result1(i) = mean(test);
end
end
plot(result1, 'b'),axis(v)
%result
%plot(result, 'r-',result, 'b:');
%axis([0,10000,-0.6,0.6])