开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 行业应用 > 邮电通讯系统 > 查看源码
Rayleigh_Doppler_multiPath.m
资源名称：多径Rayleigh信道.rar　[点击查看]
上传用户：qin0369
上传日期：2013-03-07
资源大小：9k
文件大小：8k
源码类别：

邮电通讯系统

开发平台：

Matlab

Rayleigh_Doppler_multiPath.m：源码内容
  1. function y = Rayleigh_Doppler_multiPath(fc,v,startT,endT,deltaT,fchip,delayTime,averagePower)
  2. %He jian, 2005.3
  3. %产生multipath Rayleigh分布(Doppler Shift),基于Clarke模型
  4. %return 复变量
  6. %fc=2000;%载频(MHz)
  7. %v=50;%绝对时速(km/h)
  8. % startT,endT(s)：分别表示信道仿真的开始时间、终止时间，通常startT=0，endT=1s，
  9. % deltaT(ms)：时间间隔，通常deltaT=1ms
  10. %fchip=1.28;%chip速率(Mchip/s)
  11. %delayTime=[0,781,1563,2344];%ns(10^-9s)
  12. %averagePower=[0,-3,-6,-9];%dB
  14. %method_flag = 3; %1,按ns delay，运算量大(目前不考虑该方法！)
  15.                   %2,按chip delay，运算量较大，在path delay不是chip时间整数倍时有误差
  16.                   %3,直接在抽样上delay，运算量小
  17. if(fchip~=0)
  18. method_flag = 2;
  19. else
  20. method_flag = 3;
  21. end
  23. if (method_flag==1)
  24. method_flag_str = '第1种方法，按ns delay';
  25. elseif (method_flag==2)
  26. method_flag_str = '第2种方法，按chip delay';
  27. elseif (method_flag==3)
  28. method_flag_str = '第3种方法，在抽样上delay';
  29. else
  30. end
  32. %tic;%Start timer of the simulation
  34. averagePower=10.^(averagePower/10);
  35. averagePower=averagePower/sum(averagePower);%将信道做归一化！
  37. %get number of all paths
  38. Np=length(delayTime);%delayTime与averagePower维数必须一致
  40. if (method_flag==1)%对Np条径进行delay,1st method:按ns delay，运算量大
  41. elseif (method_flag==2)%对Np条径进行delay,2nd method:按chip delay，运算量小，在path delay不是chip时间整数倍时有误差
  42. Tchip=1/(fchip*10^6);%一个chip占用时间(s)
  43. %将延时折算成延chip!!
  44. delayChip=round(delayTime*10^-9/Tchip)
  45. %get Np rayleigh path,deltaT按一个chip周期
  46. deltaT=Tchip*1000;%ms
  47. r0 = zeros(Np,length([startT:deltaT/1000:endT]));
  48. for n=1:1:Np
  49. ['==========',method_flag_str,'==========',num2str(Np),'径--- 第',num2str(n),'径==========']
  51. %tic;
  52. tempr = sqrt(averagePower(n))*Rayleigh_Doppler_singlePath(fc,v,startT,endT,deltaT);
  53. %tempr = Rayleigh_Doppler_singlePath(fc,v,startT,endT,deltaT);
  55. %delay chip...
  56. r0(n,:) = [zeros(delayChip(n),1);tempr(1:length(tempr)-delayChip(n))]';
  57. %r0(n,:)=abs(r0(n,:)).^2 * averagePower(n);%Np条径
  58. clear tempr;
  60. %disp(['one rayleigh channel time: ' num2str(toc) '秒']);
  61. end
  62. rm = sum(r0)';
  63. clear r0;
  64. elseif (method_flag==3)
  65. %此时，delayTime以最晚一个径为基准，条件参数中是以第一个径为准的！
  66. for n=1:1:Np
  67. delayTime(n) = delayTime(Np) - delayTime(n);
  68. end
  70. for n=1:1:Np
  71. lenT(n) = length([startT+delayTime(n)*10^-9:deltaT/1000:endT]);
  72. end
  73. r0 = zeros(Np,min(lenT));
  74. for n=1:1:Np
  75. ['==========',method_flag_str,'==========',num2str(Np),'径--- 第',num2str(n),'径==========']
  77. %tic;
  78. tempr = sqrt(averagePower(n))*Rayleigh_Doppler_singlePath(fc,v,startT+delayTime(n)*10^-9,endT,deltaT);
  79. %tempr = Rayleigh_Doppler_singlePath(fc,v,startT+delayTime(n)*10^-9,endT,deltaT);
  81. r0(n,:) = tempr([1:min(lenT)])';
  82. %r0(n,:)=abs(r0(n,:)).^2 * averagePower(n);%Np条径
  83. clear tempr;
  85. %disp(['one rayleigh channel time: ' num2str(toc) '秒']);
  86. end
  87. rm = sum(r0)';
  88. clear r0;
  90. startT = startT+max(delayTime)*10^-9;
  91. ['起始时间 为 startT,',num2str(startT),' s!']
  92. else
  93. end
  95. plot_flag = 0;   %1：需要plot，0：not plot
  97. if plot_flag==1
  98. fs = 1000/deltaT; %fs = fchip*10^6;
  99. c=3*10^8;%光速(m/s)
  100. fmax = (fc*10^6)*(v*10^3/3600)/c; % Max Doppler Shift (Hz)
  102. f_lim_range = [-fmax*2,fmax*2];
  104. %功率谱估计
  105. Nfft = 2^4;
  106. while(Nfft)
  107. if (Nfft < length(rm))
  108. Nfft = 2*Nfft;
  109. else
  110. break;
  111. end
  112. end
  113. Nfft = Nfft/2;%让数据长度为2的幂，又不超出采样长度
  114. r2 = rm(1:Nfft);
  116. Power_dB = 20*log10(abs(r2));% to dB!
  118. figure;
  119. subplot(2,2,1);
  120. 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值');
  121. %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值');
  122. legend(['E(r^2)=',num2str(10*log10(sum(abs(r2).^2)/length(r2)),'%.2f'),' dB']);
  123. clear Power_dB;
  125. subplot(2,2,2);
  126. clear psd_matlab;clear f_matlab;
  127. [psd_matlab,f_matlab] = pwelch(r2,[],[],'twosided',Nfft,fs);%pwelch(x,window,noverlap,nfft,fs)
  128. psd_matlab = fftshift(psd_matlab);
  129. len = length(f_matlab);
  130. plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
  131. %doppler shift
  132. hold on;
  133. sigma_u4 = sqrt(1/2);fm4 = [-fmax*0.999:fmax/100:fmax*0.999];fc4 = 0;
  134. Sf4 = 1.5*sigma_u4/(pi*fmax).*1./(sqrt(1-((fm4-fc4)./fmax).^2));
  135. 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);
  136. legend('仿真值','单径理论值');
  137. xlim(f_lim_range);grid;
  138. title('pwelch(),Welch Method');xlabel('Hz');ylabel('dB/Hz');
  140. subplot(2,2,3);
  141. clear psd_matlab;clear f_matlab;
  142. [psd_matlab,f_matlab] = pmtm(r2,4,'twosided',Nfft,fs);
  143. psd_matlab = fftshift(psd_matlab);
  144. len = length(f_matlab);
  145. plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
  146. %doppler shift
  147. hold on;
  148. 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);
  149. legend('仿真值','单径理论值');
  150. xlim(f_lim_range);grid;
  151. title('pmtm(),Multitaper method(MTM)');xlabel('Hz');ylabel('dB/Hz');
  153. subplot(2,2,4);
  154. clear psd_matlab;clear f_matlab;
  155. [psd_matlab,f_matlab] = pyulear(r2,round(Nfft/20),'twosided',Nfft,fs);
  156. psd_matlab = fftshift(psd_matlab);
  157. len = length(f_matlab);
  158. plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
  159. %doppler shift
  160. hold on;
  161. 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);
  162. legend('仿真值','单径理论值');
  163. xlim(f_lim_range);grid;
  164. title('pyulear(),Yule-Walker AR Method');xlabel('Hz');ylabel('dB/Hz');
  166. clear r2;
  167. %  yw = abs(fftshift(fft(r2))).^2/length(r2);
  168. %  clear r2;
  169. % 
  170. %  len = length(yw);
  171. %  f_range = (-len/2:len/2-1)/len*fs; %[-fs/2:1/(endT-startT):fs/2];%(0:len-1)/len*fs;
  172. %  subplot(2,2,2);plot(f_range,10*log10(yw));grid;xlim([-fmax*1.2,fmax*1.2]);title('周期图法 功率谱');xlabel('频率(Hz)');ylabel('功率谱(dB)');
  173. %  yw2=yw;f_less=find(f_range<0);f_more=find(f_range>fmax);yw2([f_less,f_more])=[];
  174. %  %legend(['(0,fmax)积分功率=',num2str(10*log10(sum(yw2)*fs/len/fmax),'%.2f'),' dB/Hz']);
  175. %  legend(['(0,fmax)积分功率=',num2str(10*log10(mean(yw2)),'%.2f'),' dB/Hz']);
  176. %  clear yw2;
  177. % 
  178. %  subplot(2,2,3);plot(f_range,10*log10(yw));xlim([-fmax*4,fmax*4]);grid;title(['周期图法 功率谱']);
  179. %  xlabel('频率(Hz)');ylabel('功率谱(dB)');
  180. %  clear yw;
  181. % 
  182. %  %======= Welch K from 2 to 5 使频域不至于展开过宽，而分辨不清!=======
  183. %  Kmax = 3; K=Kmax+1;
  184. % 
  185. %  L = 2^4; %每段数据长度,2的幂
  186. %  if (1.5*L>=length(rm))
  187. %  'Welch: 数据总长度应> 1.5*L!'
  188. %  return;
  189. %  end
  190. % 
  191. %  while (K>Kmax)
  192. %  Lmax = floor(length(rm)*2/L)/2*L; %需要从rm中提取的数据总长度
  193. %  K = Lmax*2/L-1; %数据分段数
  194. %  if (K <= Kmax)
  195. %  break;
  196. %  else
  197. %  L = 2*L;
  198. %  end
  199. %  end
  200. % 
  201. %  w_hn = hanning(L); 
  202. %  Pw = [];
  203. %  for k=1:1:K
  204. %  Pw(k,:) = (abs(fftshift(fft(w_hn.*rm(1+(k-1)*L/2:L+(k-1)*L/2)))).^2)';
  205. %  end
  206. %  Pw = sum(Pw)/(norm(w_hn)^2*K);
  207. %  f_range = (-L/2:L/2-1)/L*fs;
  208. % 
  209. %  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)');
  210. %  f_less=find(f_range<0);f_more=find(f_range>fmax);Pw([f_less,f_more])=[];
  211. %  %legend(['(0,fmax)积分功率=',num2str(10*log10(sum(Pw)*fs/L/fmax),'%.2f'),' dB/Hz']);
  212. %  legend(['(0,fmax)积分功率=',num2str(10*log10(mean(Pw)),'%.2f'),' dB/Hz']);
  213. %  clear Pw;
  214. end
  216. y=rm(:);