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

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

邮电通讯系统

开发平台：

Matlab

Rayleigh_Doppler_singlePath.m：源码内容
  1. function y = Rayleigh_Doppler_singlePath(fc,v,startT,endT,deltaT)
  2. %He jian, 2005.3
  3. %产生单径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. if (endT-startT<=sectionTime/1000)
  11. '计算时间段必须小于总时长!'
  12. return;
  13. end
  14. if (sectionTime>0&&deltaT>=sectionTime)
  15. 'deltaT要小于时间段!'
  16. return;
  17. end
  19. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  20. %6种产生单径Rayleigh方法，参考如下：
  21. %1&2: A deterministic digital simulation model for Suzuki processes with application to a shadowed Rayleigh channel
  22. %3  : Rayleigh Channel Fading Simulator：Problems and Solutions
  23. %4  : 频域滤波方法参见rayleigh_Filter_Model.m
  24. %5  : Improved Models for the Generation of Multiple Uncorrelated Rayleigh Fading Waveforms
  25. %6  : Simulation Models With Correct Statistical Properties for Rayleigh Fading Channels
  26. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  28. method_flag = 6; %1,Suzuki-MEA;?只有正频率部分 对该方法(一次计算部分)进行了优化处理!
  29.  %2,Suzuki-MED;?好象这个有点问题，出来的pdf与rayleigh比较，波动较大，不如MEA！
  30.  %3,New model(t);对该方法(一次计算部分)进行了优化处理!
  31.  %4,Filter model(f); rayleigh_Filter_Model.m
  32.  %5,New model(t) IEEE.2002.06;对该方法(一次计算部分)进行了优化处理!
  33.  %6,New model(t) IEEE.2003.06;对该方法进行了优化处理!
  34.  
  35. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  36. %检查该方法的结果是否为rayleigh的pdf！(幅度|Z(t)|的PDF)
  37. %clear;r=Rayleigh_Doppler_singlePath(2000,30,0,20,1,0);test_rayleigh_pdf(r);
  38. %检查其相位分布？？
  39. %phase=angle(xcorr(r))*180/pi;figure;subplot(2,1,1);plot(phase);subplot(2,1,2);hist(phase,50);
  40. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  42. if (method_flag==1)
  43. method_str = ['Suzuki-MEA'];
  44. elseif (method_flag==2)
  45. method_str = ['Suzuki-MED'];
  46. elseif (method_flag==3)
  47. method_str = ['New model(t)'];
  48. elseif (method_flag==4)
  49. method_str = ['Filter model(f)'];
  50. elseif (method_flag==5)
  51. method_str = ['New model(t) IEEE.2002.06'];
  52. elseif (method_flag==6)
  53. method_str = ['New model(t) IEEE.2003.06'];
  54. else
  55. end
  57. sigma_u = sqrt(1/2);%为rayleigh分布pdf的参数
  59. c=3*10^8;%光速(m/s)
  60. fmax = (fc*10^6)*(v*10^3/3600)/c; % Max Doppler Shift (Hz)
  62. fs=1000/deltaT;%抽样频率,if fs=1000-->模拟时间间隔=1/1000秒，即1ms
  64. rand('state',sum(100*clock));
  66. %如果涉及fft等计算，在Nfm点数为2的幂次方时，计算效率最高！
  67. %确保2*fmax的频率范围内，抽样点数满足2的幂次方要求！同时又要超过抽样点数！
  68. %   |---------------|----------------------........----------------------------|
  69. %   0              2*fmax                                                     fs
  70. %   |<---- Nfm ---->|                                                          |     
  71. %   |<----------------------------- Nfs -------------------------------------->|
  73. tic;
  75. N = 32;%模拟均匀到达的rays数
  77. if(method_flag==1)%MEA
  78. fd = fmax*sin((1:N)*pi/2/N); %cos(2*pi*((1:N)/N))*fmax; 
  79. phase = unifrnd(0,2*pi,N,1);  %phase = 2*pi*((1:N)/N);%相位 0-2*pi 均匀分布 
  80. %phase = unifrnd(-pi,pi,N,1);   
  82. Nt = length([startT:deltaT/1000:endT]);
  83. %temp = ones(1,Nt);
  84. %phase = phase(:)*temp;
  86. step = 0;
  87. r1 = zeros(length(startT:deltaT/1000:endT),1);
  88. for t = startT:deltaT/1000:endT
  89. step = step + 1;
  90. r1(step) = sum((exp(i*(2*pi*fd(:)*t + phase(:)))))*sigma_u*sqrt(2/N);%幅度
  91. end
  92. %t = [startT:deltaT/1000:endT];
  93. %r1 = sum((exp(i*(2*pi*fd(:)*t + phase))))*sigma_u*sqrt(2/N);%幅度
  94. %r1 = r1(:);
  95. elseif(method_flag==2)%MED
  96. fd = fmax*(2*(1:N)-1)/(2*N);%fd = fmax*(-N/2:N/2-1)/N; 
  97. phase = unifrnd(0,2*pi,N,1);  %phase = 2*pi*((1:N)/N);%相位 0-2*pi 均匀分布 
  98. cn = sqrt(4*sigma_u*sigma_u*(asin((1:N)/N)-asin(((1:N)-1)/N))/pi);
  100. Nt = length([startT:deltaT/1000:endT]);
  102. step = 0;
  103. for t = startT:deltaT/1000:endT
  104. step = step + 1;
  105. r1(step) = sum(exp(i*(2*pi*fd(:)*t + phase(:))).*cn(:));%幅度
  106. end
  107. r1 = r1(:);
  108. elseif(method_flag==3)%New model(t)
  109. M3 = 16; %>=8即可！
  110. N3 = 4*M3 + 2;
  112. theta = unifrnd(-pi,pi,M3+1,1);  %[-pi,pi)均匀分布
  113. fi = unifrnd(-pi,pi,M3+1,1);  %[-pi,pi)均匀分布
  114. ys = unifrnd(-pi,pi,M3+1,1);  %[-pi,pi)均匀分布
  116. wd = 2*pi*fmax;
  117. wn = wd*cos(2*pi*[0:M3]'/N3+theta(:)/N3);
  119. step = 0;
  120. r1 = zeros(length(startT:deltaT/1000:endT),1);
  121. for t = startT:deltaT/1000:endT
  122. step = step + 1;
  123. Xc = 2/sqrt(N3)*((sqrt(2)*cos(ys(1))*cos(wn(1)*t+fi(1)))+sum(2*cos(ys(2:end)).*cos(wn(2:end)*t+fi(2:end))));
  124. Xs = 2/sqrt(N3)*((sqrt(2)*sin(ys(1))*cos(wn(1)*t+fi(1)))+sum(2*sin(ys(2:end)).*cos(wn(2:end)*t+fi(2:end))));
  125. r1(step) = Xc + i*Xs;%幅度
  126. r1(step) = r1(step) * sigma_u;
  127. end
  128. %r1 = r1(:);
  129. elseif(method_flag==4)%Filter model(f)
  130. r1 = rayleigh_Filter_Model(fmax,fs,(endT-startT)*fs+1);
  131. r1 = r1(:);
  132. elseif(method_flag==5)%New model(t) IEEE.2002.06
  133. M5 = 32;%>=8
  135. %theta = unifrnd(-pi,pi,M5,1);  %[-pi,pi)均匀分布
  136. theta = unifrnd(-pi,pi,1,1);  %[-pi,pi)均匀分布
  138. fi = unifrnd(-pi,pi,M5,1);  %[-pi,pi)均匀分布
  139. ys = unifrnd(-pi,pi,M5,1);  %[-pi,pi)均匀分布
  140. a5 = (2*pi*[1:M5]'-pi+theta)/(4*M5);
  141. step = 0;
  142. r1 = zeros(length(startT:deltaT/1000:endT),1);
  143. for t = startT:deltaT/1000:endT
  144. step = step + 1;
  145. Zc = sqrt(2/M5)* sum(cos(2*pi*fmax*t*cos(a5)+fi));
  146. Zs = sqrt(2/M5)* sum(cos(2*pi*fmax*t*sin(a5)+ys));
  147. r1(step) = Zc + i*Zs;%幅度
  148. r1(step) = r1(step) * sigma_u;
  149. end
  150. %r1 = r1(:);
  151. elseif(method_flag==6)%New model(t) IEEE.2003.06
  152. '一次计算中... ...'
  153. M6 = 32;%>=8
  154. %theta = unifrnd(-pi,pi,M6,1);  %[-pi,pi)均匀分布
  155. %fi = unifrnd(-pi,pi,M6,1);  %[-pi,pi)均匀分布
  156. theta = unifrnd(-pi,pi,1,1);  %[-pi,pi)均匀分布
  157. fi = unifrnd(-pi,pi,1,1);  %[-pi,pi)均匀分布
  159. ys = unifrnd(-pi,pi,M6,1);  %[-pi,pi)均匀分布
  161. a6 = (2*pi*[1:M6]'-pi+theta)/(4*M6);
  162. step = 0;
  164. r1 = zeros(length(startT:deltaT/1000:endT),1);
  165. for t = startT:deltaT/1000:endT
  166. step = step + 1;
  167. Xc = sqrt(4/M6)* sum(cos(ys).*cos(2*pi*fmax*t*cos(a6)+fi));
  168. Xs = sqrt(4/M6)* sum(sin(ys).*cos(2*pi*fmax*t*cos(a6)+fi));
  169. r1(step) = Xc + i*Xs;%幅度
  170. r1(step) = r1(step) * sigma_u;
  171. end
  172. %r1 = r1(:);
  173. else
  174. end
  175. disp(['one rayleigh channel time: ' num2str(toc) '秒']);
  177. plot_flag = 0;   %2：-包括own的Welch法,周期图法(对高采样率优化plot);
  178.                  %       matlab的periodogram,pwelch,pmtm,pyulear;
  179.  %   -还可以计算积分功率;
  180.                  %1：需要一般plot(只有周期图法)
  181.                  %0：not plot
  182.                  
  183. if(method_flag==1)
  184. f_lim_range = [-fmax*0.2,fmax*2];
  185. elseif(method_flag==2)
  186. f_lim_range = [-fmax*0.2,fmax*2];
  187. elseif(method_flag==3)
  188. f_lim_range = [-fmax*2,fmax*2];
  189. elseif(method_flag==4)
  190. f_lim_range = [-fmax*2,fmax*2];
  191. elseif(method_flag==5)
  192. f_lim_range = [-fmax*2,fmax*2];
  193. elseif(method_flag==6)
  194. f_lim_range = [-fmax*2,fmax*2];
  195. end
  197. if plot_flag==1
  198. Power_dB = 10*log10(abs(r1).^2);
  199. subplot(2,2,1);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,',method_str]);xlabel('ms');ylabel('dB值');
  201. %功率谱估计
  202. yw = fft(r1);
  203. yw=fftshift(yw);
  205. yw = abs(yw).^2/length(yw);
  206. len = length(yw);
  207. f_range = [-fs/2:1/(endT-startT):fs/2];%(0:len-1)/len*fs;
  208. subplot(2,2,2);plot(f_range,10*log10(yw));grid;xlim(f_lim_range);title('周期图法 功率谱 abs()^2/N');xlabel('频率(Hz)');ylabel('功率谱(dB)');
  210. rt = xcorr(r1,r1);
  211. subplot(2,2,3);plot(10*log10(abs(rt)/length(rt)));grid;axis tight;title(['自相关性']);
  213. clear yw;
  214. yw = xcorr(r1,r1);
  215. len = length(yw);
  216. yw = fft(yw);yw=fftshift(yw);
  218. Pyw = abs(yw)/length(yw);%.^2/length(yw);
  219. subplot(2,2,4);
  220. plot((-len/2:len/2-1)/len*fs,10*log10(Pyw));grid;xlim(f_lim_range);title(['自相关法 功率谱,plot_flag=',num2str(plot_flag)]);xlabel('频率(Hz)');ylabel('功率谱(dB)');
  221. elseif plot_flag==2
  222. %功率谱估计
  223. fs = 1000/deltaT; 
  224. Nfft = 2^4;
  226. while(Nfft)
  227. if (Nfft < length(r1))
  228. Nfft = 2*Nfft;
  229. else
  230. break;
  231. end
  232. end
  234. Nfft = Nfft/2;%让数据长度为2的幂，又不超出采样长度
  235. r2 = r1(1:Nfft);
  237. Power_dB = 20*log10(abs(r2));%dB
  238. figure;
  239. 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,',method_str]);xlabel('ms');ylabel('dB值');
  240. legend(['E(r^2)=',num2str(10*log10(sum(abs(r2).^2)/length(r2)),'%.2f'),' dB']);
  241. clear Power_dB;
  244. %  yw = abs(fftshift(fft(r2))).^2/length(r2);
  245. % 
  246. %  len = length(yw);
  247. %  f_range = (-len/2:len/2-1)/len*fs; %[-fs/2:1/(endT-startT):fs/2];%(0:len-1)/len*fs;
  248. %  subplot(2,2,2);plot(f_range,10*log10(yw));grid;xlim(f_lim_range);title('周期图法 功率谱');xlabel('频率(Hz)');ylabel('功率谱(dB)');
  249. %  yw2=yw;f_less=find(f_range<0);f_more=find(f_range>fmax);yw2([f_less,f_more])=[];
  250. %  %size(yw2),fmax*len/fs
  251. %  %legend(['(0,fmax)mean=',num2str(10*log10(sum(yw2)*fs/len/fmax),'%.2f'),' dB/Hz']);
  252. %  legend(['(0,fmax)mean=',num2str(10*log10(mean(yw2)),'%.2f'),' dB/Hz']);
  253. %  clear yw2;
  254. % 
  255. %  subplot(2,2,3);plot(f_range,10*log10(yw));grid;xlim([-fmax*4,fmax*4]);title(['周期图法 功率谱']);xlabel('频率(Hz)');ylabel('功率谱(dB)');
  256. %  clear yw;
  257. % 
  258. %  %======= Welch K from 2 to 5 使频域不至于展开过宽，而分辨不清!=======
  259. %  Kmax = 3; K=Kmax+1;
  260. % 
  261. %  L = 2^4; %每段数据长度,2的幂
  262. %  if (1.5*L>=length(r2))
  263. %  'Welch: 数据总长度应> 1.5*L!'
  264. %  return;
  265. %  end
  266. % 
  267. %  while (K>Kmax)
  268. %  Lmax = floor(length(r2)*2/L)/2*L; %需要从r2中提取的数据总长度
  269. %  K = Lmax*2/L-1; %数据分段数
  270. %  if (K <= Kmax)
  271. %  if (K==1)
  272. %  K = 2;
  273. %  L = L/2;
  274. %  end
  275. %  break;
  276. %  else
  277. %  L = 2*L;
  278. %  end
  279. %  end
  280. % 
  281. %  w_hn = hanning(L); 
  282. %  Pw = [];
  283. %  for k=1:1:K
  284. %  Pw(k,:) = (abs(fftshift(fft(w_hn.*r2(1+(k-1)*L/2:L+(k-1)*L/2)))).^2)';
  285. %  end
  286. % 
  287. %  Pw = sum(Pw)/(norm(w_hn)^2*K);
  288. %  f_range = (-L/2:L/2-1)/L*fs;
  289. % 
  290. %  subplot(2,2,4);plot(f_range,10*log10(Pw));grid;xlim(f_lim_range);title(['Welch法 功率谱,K=',num2str(K),',L=2^',num2str(log2(L)),',plot_flag=',num2str(plot_flag)]);xlabel('频率(Hz)');ylabel('功率谱(dB)');
  291. %  f_less=find(f_range<0);f_more=find(f_range>fmax);Pw([f_less,f_more])=[];
  292. %  %legend(['(0,fmax)mean=',num2str(10*log10(sum(Pw)*fs/L/fmax),'%.2f'),' dB/Hz']);
  293. %  legend(['(0,fmax)mean=',num2str(10*log10(mean(Pw)),'%.2f'),' dB/Hz']);
  294. %  clear Pw;
  295. % 
  296. %  figure;
  297. %  subplot(2,2,1);
  298. %  [psd_matlab,f_matlab] = periodogram(r2,[],'twosided',Nfft,fs);%
  299. %  psd_matlab = fftshift(psd_matlab);
  300. %  len = length(f_matlab);
  301. %  plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));xlim(f_lim_range);grid;
  302. %  title(['periodogram(),',num2str(fc), 'MHz,',num2str(v), 'km/h,Max Doppler=',num2str(fmax,'%.2f'),'Hz,']);
  303. %  xlabel('Hz');ylabel('dB/Hz');
  305. subplot(2,2,2);
  306. clear psd_matlab;clear f_matlab;
  307. [psd_matlab,f_matlab] = pwelch(r2,[],[],'twosided',Nfft,fs);%pwelch(x,window,noverlap,nfft,fs)
  308. psd_matlab = fftshift(psd_matlab);
  309. len = length(f_matlab);
  310. plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
  312. %doppler shift
  313. hold on;
  314. sigma_u4 = sqrt(1/2);fm4 = [-fmax*0.999:fmax/100:fmax*0.999];fc4 = 0;
  315. Sf4 = 1.5*sigma_u4/(pi*fmax).*1./(sqrt(1-((fm4-fc4)./fmax).^2));
  316. 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);
  317. legend('仿真值','理论值');
  319. xlim(f_lim_range);grid;
  320. title('pwelch(),Welch Method');xlabel('Hz');ylabel('dB/Hz');
  322. subplot(2,2,3);
  323. clear psd_matlab;clear f_matlab;
  324. [psd_matlab,f_matlab] = pmtm(r2,4,'twosided',Nfft,fs);
  325. psd_matlab = fftshift(psd_matlab);
  326. len = length(f_matlab);
  327. plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
  329. %doppler shift
  330. hold on;
  331. 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);
  332. legend('仿真值','理论值');
  334. xlim(f_lim_range);grid;
  335. title('pmtm(),Multitaper method(MTM)');xlabel('Hz');ylabel('dB/Hz');
  337. subplot(2,2,4);
  338. clear psd_matlab;clear f_matlab;
  339. [psd_matlab,f_matlab] = pyulear(r2,round(Nfft/20),'twosided',Nfft,fs);
  340. psd_matlab = fftshift(psd_matlab);
  341. len = length(f_matlab);
  342. plot([-flipud(f_matlab(2:len/2+1));f_matlab(1:len/2)],10*log10(psd_matlab));
  344. %doppler shift
  345. hold on;
  346. 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);
  347. legend('仿真值','理论值');
  349. xlim(f_lim_range);grid;
  350. title('pyulear(),Yule-Walker AR Method');xlabel('Hz');ylabel('dB/Hz');
  352. clear r2;
  353. else
  354. end
  356. y=r1(:);%信道复变量