verticr_recursrve_snapshot.zip_snapshot_多经信道_无线信道

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