rayleigh_doppler_multipath.m

服从瑞丽分布的多普勒频移的MATLAB仿真

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(:);