function y = Rayleigh_Doppler_multiPath(fc,v,startT,endT,deltaT,fchip,delayTime,averagePower)
%He jian, 2005.3
%产生multipath Rayleigh分布(Doppler Shift),基于Clarke模型
%return 复变量
%fc=2000;%载频(MHz)
%v=50;%绝对时速(km/h)
% startT,endT(s)：分别表示信道仿真的开始时间、终止时间，通常startT=0，endT=1s，
% deltaT(ms)：时间间隔，通常deltaT=1ms
%fchip=1.28;%chip速率(Mchip/s)
%delayTime=[0,781,1563,2344];%ns(10^-9s)
%averagePower=[0,-3,-6,-9];%dB
%method_flag = 3; %1,按ns delay，运算量大(目前不考虑该方法！)
%2,按chip delay，运算量较大，在path delay不是chip时间整数倍时有误差
%3,直接在抽样上delay，运算量小
if(fchip~=0)
method_flag = 2;
else
method_flag = 3;
end
if (method_flag==1)
method_flag_str = '第1种方法，按ns delay';
elseif (method_flag==2)
method_flag_str = '第2种方法，按chip delay';
elseif (method_flag==3)
method_flag_str = '第3种方法，在抽样上delay';
else
end
%tic;%Start timer of the simulation
averagePower=10.^(averagePower/10);
averagePower=averagePower/sum(averagePower);%将信道做归一化！
%get number of all paths
Np=length(delayTime);%delayTime与averagePower维数必须一致
if (method_flag==1)%对Np条径进行delay,1st method:按ns delay，运算量大
elseif (method_flag==2)%对Np条径进行delay,2nd method:按chip delay，运算量小，在path delay不是chip时间整数倍时有误差
Tchip=1/(fchip*10^6);%一个chip占用时间(s)
%将延时折算成延chip!!
delayChip=round(delayTime*10^-9/Tchip)
%get Np rayleigh path,deltaT按一个chip周期
deltaT=Tchip*1000;%ms
r0 = zeros(Np,length([startT:deltaT/1000:endT]));
for n=1:1:Np
['==========',method_flag_str,'==========',num2str(Np),'径--- 第',num2str(n),'径==========']
%tic;
tempr = sqrt(averagePower(n))*Rayleigh_Doppler_singlePath(fc,v,startT,endT,deltaT);
%tempr = Rayleigh_Doppler_singlePath(fc,v,startT,endT,deltaT);
%delay chip...
r0(n,:) = [zeros(delayChip(n),1);tempr(1:length(tempr)-delayChip(n))]';
%r0(n,:)=abs(r0(n,:)).^2 * averagePower(n);%Np条径
clear tempr;
%disp(['one rayleigh channel time: ' num2str(toc) '秒']);
end
rm = sum(r0)';
clear r0;
elseif (method_flag==3)
%此时，delayTime以最晚一个径为基准，条件参数中是以第一个径为准的！
for n=1:1:Np
delayTime(n) = delayTime(Np) - delayTime(n);
end
for n=1:1:Np
lenT(n) = length([startT+delayTime(n)*10^-9:deltaT/1000:endT]);
end
r0 = zeros(Np,min(lenT));
for n=1:1:Np
['==========',method_flag_str,'==========',num2str(Np),'径--- 第',num2str(n),'径==========']
%tic;
tempr = sqrt(averagePower(n))*Rayleigh_Doppler_singlePath(fc,v,startT+delayTime(n)*10^-9,endT,deltaT);
%tempr = Rayleigh_Doppler_singlePath(fc,v,startT+delayTime(n)*10^-9,endT,deltaT);
r0(n,:) = tempr([1:min(lenT)])';
%r0(n,:)=abs(r0(n,:)).^2 * averagePower(n);%Np条径
clear tempr;
%disp(['one rayleigh channel time: ' num2str(toc) '秒']);
end
rm = sum(r0)';
clear r0;
startT = startT+max(delayTime)*10^-9;
['起始时间 为 startT,',num2str(startT),' s!']
else
end
plot_flag = 0; %1：需要plot，0：not plot
if plot_flag==1
fs = 1000/deltaT; %fs = fchip*10^6;
c=3*10^8;%光速(m/s)
fmax = (fc*10^6)*(v*10^3/3600)/c; % Max Doppler Shift (Hz)
f_lim_range = [-fmax*2,fmax*2];
%功率谱估计
Nfft = 2^4;
while(Nfft)
if (Nfft < length(rm))
Nfft = 2*Nfft;
else
break;
end
end
Nfft = Nfft/2;%让数据长度为2的幂，又不超出采样长度
r2 = rm(1:Nfft);
Power_dB = 20*log10(abs(r2));% to dB!
figure;
subplot(2,2,1);
plot([startT*1000:deltaT:startT*1000+deltaT*(Nfft-1)],Power_dB);grid;axis tight;title([num2str(fc), 'MHz,',num2str(v), 'km/h,Max Doppler=',num2str(fmax,'%.2f'),'Hz,',num2str(Np),'条径']);xlabel('ms');ylabel('dB值');
%plot([startT*1000:deltaT:endT*1000],Power_dB);grid;axis tight;title([num2str(fc), 'MHz,',num2str(v), 'km/h,Max Doppler=',num2str(fmax,'%.2f'),'Hz,',num2str(Np),'条径']);xlabel('ms');ylabel('dB值');
legend(['E(r^2)=',num2str(10*log10(sum(abs(r2).^2)/length(r2)),'%.2f'),' dB']);
clear Power_dB;
subplot(2,2,2);
clear psd_matlab;clear f_matlab;
[psd_matlab,f_matlab] = pwelch(r2,[],[],'twosided',Nfft,fs);%pwelch(x,window,noverlap,nfft,fs)
psd_matlab = fftshift(psd_matlab);
len = length(f_matlab);
plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
%doppler shift
hold on;
sigma_u4 = sqrt(1/2);fm4 = [-fmax*0.999:fmax/100:fmax*0.999];fc4 = 0;
Sf4 = 1.5*sigma_u4/(pi*fmax).*1./(sqrt(1-((fm4-fc4)./fmax).^2));
plot(fm4,10*log10(Sf4),'-.r',min(fm4).*ones(1,2),[min(10*log10(psd_matlab)),max(10*log10(Sf4))],'-.r',max(fm4).*ones(1,2),[min(10*log10(psd_matlab)),max(10*log10(Sf4))],'-.r','LineWidth',1.5);
legend('仿真值','单径理论值');
xlim(f_lim_range);grid;
title('pwelch(),Welch Method');xlabel('Hz');ylabel('dB/Hz');
subplot(2,2,3);
clear psd_matlab;clear f_matlab;
[psd_matlab,f_matlab] = pmtm(r2,4,'twosided',Nfft,fs);
psd_matlab = fftshift(psd_matlab);
len = length(f_matlab);
plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
%doppler shift
hold on;
plot(fm4,10*log10(Sf4),'-.r',min(fm4).*ones(1,2),[min(10*log10(psd_matlab)),max(10*log10(Sf4))],'-.r',max(fm4).*ones(1,2),[min(10*log10(psd_matlab)),max(10*log10(Sf4))],'-.r','LineWidth',1.5);
legend('仿真值','单径理论值');
xlim(f_lim_range);grid;
title('pmtm(),Multitaper method(MTM)');xlabel('Hz');ylabel('dB/Hz');
subplot(2,2,4);
clear psd_matlab;clear f_matlab;
[psd_matlab,f_matlab] = pyulear(r2,round(Nfft/20),'twosided',Nfft,fs);
psd_matlab = fftshift(psd_matlab);
len = length(f_matlab);
plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
%doppler shift
hold on;
plot(fm4,10*log10(Sf4),'-.r',min(fm4).*ones(1,2),[min(10*log10(psd_matlab)),max(10*log10(Sf4))],'-.r',max(fm4).*ones(1,2),[min(10*log10(psd_matlab)),max(10*log10(Sf4))],'-.r','LineWidth',1.5);
legend('仿真值','单径理论值');
xlim(f_lim_range);grid;
title('pyulear(),Yule-Walker AR Method');xlabel('Hz');ylabel('dB/Hz');
clear r2;
% yw = abs(fftshift(fft(r2))).^2/length(r2);
% clear r2;
%
% len = length(yw);
% f_range = (-len/2:len/2-1)/len*fs; %[-fs/2:1/(endT-startT):fs/2];%(0:len-1)/len*fs;
% subplot(2,2,2);plot(f_range,10*log10(yw));grid;xlim([-fmax*1.2,fmax*1.2]);title('周期图法 功率谱');xlabel('频率(Hz)');ylabel('功率谱(dB)');
% yw2=yw;f_less=find(f_range<0);f_more=find(f_range>fmax);yw2([f_less,f_more])=[];
% %legend(['(0,fmax)积分功率=',num2str(10*log10(sum(yw2)*fs/len/fmax),'%.2f'),' dB/Hz']);
% legend(['(0,fmax)积分功率=',num2str(10*log10(mean(yw2)),'%.2f'),' dB/Hz']);
% clear yw2;
%
% subplot(2,2,3);plot(f_range,10*log10(yw));xlim([-fmax*4,fmax*4]);grid;title(['周期图法 功率谱']);
% xlabel('频率(Hz)');ylabel('功率谱(dB)');
% clear yw;
%
% %======= Welch K from 2 to 5 使频域不至于展开过宽，而分辨不清!=======
% Kmax = 3; K=Kmax+1;
%
% L = 2^4; %每段数据长度,2的幂
% if (1.5*L>=length(rm))
% 'Welch: 数据总长度应> 1.5*L!'
% return;
% end
%
% while (K>Kmax)
% Lmax = floor(length(rm)*2/L)/2*L; %需要从rm中提取的数据总长度
% K = Lmax*2/L-1; %数据分段数
% if (K <= Kmax)
% break;
% else
% L = 2*L;
% end
% end
%
% w_hn = hanning(L);
% Pw = [];
% for k=1:1:K
% Pw(k,:) = (abs(fftshift(fft(w_hn.*rm(1+(k-1)*L/2:L+(k-1)*L/2)))).^2)';
% end
% Pw = sum(Pw)/(norm(w_hn)^2*K);
% f_range = (-L/2:L/2-1)/L*fs;
%
% subplot(2,2,4);plot(f_range,10*log10(Pw));grid;title(['Welch法 功率谱,K=',num2str(K),',L=2\^',num2str(log2(L))]);xlim([-fmax*1.2,fmax*1.2]);xlabel('频率(Hz)');ylabel('功率谱(dB)');
% f_less=find(f_range<0);f_more=find(f_range>fmax);Pw([f_less,f_more])=[];
% %legend(['(0,fmax)积分功率=',num2str(10*log10(sum(Pw)*fs/L/fmax),'%.2f'),' dB/Hz']);
% legend(['(0,fmax)积分功率=',num2str(10*log10(mean(Pw)),'%.2f'),' dB/Hz']);
% clear Pw;
end
y=rm(:);
