📄 rayleigh_doppler_singlepath.m
字号:
function y = Rayleigh_Doppler_singlePath(fc,v,startT,endT,deltaT)
%He jian, 2005.3
%产生单径Rayleigh分布(Doppler Shift),基于Clarke模型
%return 信道复变量
%fc=2000;%载频(MHz)
%v=50;%绝对时速(km/h)
% startT,endT(s)：分别表示信道仿真的开始时间、终止时间，通常startT=0，endT=1s，
% deltaT(ms)：时间间隔，通常deltaT=1ms
if (endT-startT<=sectionTime/1000)
'计算时间段必须小于总时长!'
return;
end
if (sectionTime>0&&deltaT>=sectionTime)
'deltaT要小于时间段!'
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%6种产生单径Rayleigh方法，参考如下：
%1&2: A deterministic digital simulation model for Suzuki processes with application to a shadowed Rayleigh channel
%3 : Rayleigh Channel Fading Simulator：Problems and Solutions
%4 : 频域滤波方法参见rayleigh_Filter_Model.m
%5 : Improved Models for the Generation of Multiple Uncorrelated Rayleigh Fading Waveforms
%6 : Simulation Models With Correct Statistical Properties for Rayleigh Fading Channels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
method_flag = 6; %1,Suzuki-MEA;?只有正频率部分 对该方法(一次计算部分)进行了优化处理!
%2,Suzuki-MED;?好象这个有点问题，出来的pdf与rayleigh比较，波动较大，不如MEA！
%3,New model(t);对该方法(一次计算部分)进行了优化处理!
%4,Filter model(f); rayleigh_Filter_Model.m
%5,New model(t) IEEE.2002.06;对该方法(一次计算部分)进行了优化处理!
%6,New model(t) IEEE.2003.06;对该方法进行了优化处理!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%检查该方法的结果是否为rayleigh的pdf！(幅度|Z(t)|的PDF)
%clear;r=Rayleigh_Doppler_singlePath(2000,30,0,20,1,0);test_rayleigh_pdf(r);
%检查其相位分布？？
%phase=angle(xcorr(r))*180/pi;figure;subplot(2,1,1);plot(phase);subplot(2,1,2);hist(phase,50);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (method_flag==1)
method_str = ['Suzuki-MEA'];
elseif (method_flag==2)
method_str = ['Suzuki-MED'];
elseif (method_flag==3)
method_str = ['New model(t)'];
elseif (method_flag==4)
method_str = ['Filter model(f)'];
elseif (method_flag==5)
method_str = ['New model(t) IEEE.2002.06'];
elseif (method_flag==6)
method_str = ['New model(t) IEEE.2003.06'];
else
end
sigma_u = sqrt(1/2);%为rayleigh分布pdf的参数
c=3*10^8;%光速(m/s)
fmax = (fc*10^6)*(v*10^3/3600)/c; % Max Doppler Shift (Hz)
fs=1000/deltaT;%抽样频率,if fs=1000-->模拟时间间隔=1/1000秒，即1ms
rand('state',sum(100*clock));
%如果涉及fft等计算，在Nfm点数为2的幂次方时，计算效率最高！
%确保2*fmax的频率范围内，抽样点数满足2的幂次方要求！同时又要超过抽样点数！
% |---------------|----------------------........----------------------------|
% 0 2*fmax fs
% |<---- Nfm ---->| |
% |<----------------------------- Nfs -------------------------------------->|
tic;
N = 32;%模拟均匀到达的rays数
if(method_flag==1)%MEA
fd = fmax*sin((1:N)*pi/2/N); %cos(2*pi*((1:N)/N))*fmax;
phase = unifrnd(0,2*pi,N,1); %phase = 2*pi*((1:N)/N);%相位 0-2*pi 均匀分布
%phase = unifrnd(-pi,pi,N,1);
Nt = length([startT:deltaT/1000:endT]);
%temp = ones(1,Nt);
%phase = phase(:)*temp;
step = 0;
r1 = zeros(length(startT:deltaT/1000:endT),1);
for t = startT:deltaT/1000:endT
step = step + 1;
r1(step) = sum((exp(i*(2*pi*fd(:)*t + phase(:)))))*sigma_u*sqrt(2/N);%幅度
end
%t = [startT:deltaT/1000:endT];
%r1 = sum((exp(i*(2*pi*fd(:)*t + phase))))*sigma_u*sqrt(2/N);%幅度
%r1 = r1(:);
elseif(method_flag==2)%MED
fd = fmax*(2*(1:N)-1)/(2*N);%fd = fmax*(-N/2:N/2-1)/N;
phase = unifrnd(0,2*pi,N,1); %phase = 2*pi*((1:N)/N);%相位 0-2*pi 均匀分布
cn = sqrt(4*sigma_u*sigma_u*(asin((1:N)/N)-asin(((1:N)-1)/N))/pi);
Nt = length([startT:deltaT/1000:endT]);
step = 0;
for t = startT:deltaT/1000:endT
step = step + 1;
r1(step) = sum(exp(i*(2*pi*fd(:)*t + phase(:))).*cn(:));%幅度
end
r1 = r1(:);
elseif(method_flag==3)%New model(t)
M3 = 16; %>=8即可！
N3 = 4*M3 + 2;
theta = unifrnd(-pi,pi,M3+1,1); %[-pi,pi)均匀分布
fi = unifrnd(-pi,pi,M3+1,1); %[-pi,pi)均匀分布
ys = unifrnd(-pi,pi,M3+1,1); %[-pi,pi)均匀分布
wd = 2*pi*fmax;
wn = wd*cos(2*pi*[0:M3]'/N3+theta(:)/N3);
step = 0;
r1 = zeros(length(startT:deltaT/1000:endT),1);
for t = startT:deltaT/1000:endT
step = step + 1;
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))));
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))));
r1(step) = Xc + i*Xs;%幅度
r1(step) = r1(step) * sigma_u;
end
%r1 = r1(:);
elseif(method_flag==4)%Filter model(f)
r1 = rayleigh_Filter_Model(fmax,fs,(endT-startT)*fs+1);
r1 = r1(:);
elseif(method_flag==5)%New model(t) IEEE.2002.06
M5 = 32;%>=8
%theta = unifrnd(-pi,pi,M5,1); %[-pi,pi)均匀分布
theta = unifrnd(-pi,pi,1,1); %[-pi,pi)均匀分布
fi = unifrnd(-pi,pi,M5,1); %[-pi,pi)均匀分布
ys = unifrnd(-pi,pi,M5,1); %[-pi,pi)均匀分布
a5 = (2*pi*[1:M5]'-pi+theta)/(4*M5);
step = 0;
r1 = zeros(length(startT:deltaT/1000:endT),1);
for t = startT:deltaT/1000:endT
step = step + 1;
Zc = sqrt(2/M5)* sum(cos(2*pi*fmax*t*cos(a5)+fi));
Zs = sqrt(2/M5)* sum(cos(2*pi*fmax*t*sin(a5)+ys));
r1(step) = Zc + i*Zs;%幅度
r1(step) = r1(step) * sigma_u;
end
%r1 = r1(:);
elseif(method_flag==6)%New model(t) IEEE.2003.06
'一次计算中... ...'
M6 = 32;%>=8
%theta = unifrnd(-pi,pi,M6,1); %[-pi,pi)均匀分布
%fi = unifrnd(-pi,pi,M6,1); %[-pi,pi)均匀分布
theta = unifrnd(-pi,pi,1,1); %[-pi,pi)均匀分布
fi = unifrnd(-pi,pi,1,1); %[-pi,pi)均匀分布
ys = unifrnd(-pi,pi,M6,1); %[-pi,pi)均匀分布
a6 = (2*pi*[1:M6]'-pi+theta)/(4*M6);
step = 0;
r1 = zeros(length(startT:deltaT/1000:endT),1);
for t = startT:deltaT/1000:endT
step = step + 1;
Xc = sqrt(4/M6)* sum(cos(ys).*cos(2*pi*fmax*t*cos(a6)+fi));
Xs = sqrt(4/M6)* sum(sin(ys).*cos(2*pi*fmax*t*cos(a6)+fi));
r1(step) = Xc + i*Xs;%幅度
r1(step) = r1(step) * sigma_u;
end
%r1 = r1(:);
else
end
disp(['one rayleigh channel time: ' num2str(toc) '秒']);
plot_flag = 0; %2：-包括own的Welch法,周期图法(对高采样率优化plot);
% matlab的periodogram,pwelch,pmtm,pyulear;
% -还可以计算积分功率;
%1：需要一般plot(只有周期图法)
%0：not plot
if(method_flag==1)
f_lim_range = [-fmax*0.2,fmax*2];
elseif(method_flag==2)
f_lim_range = [-fmax*0.2,fmax*2];
elseif(method_flag==3)
f_lim_range = [-fmax*2,fmax*2];
elseif(method_flag==4)
f_lim_range = [-fmax*2,fmax*2];
elseif(method_flag==5)
f_lim_range = [-fmax*2,fmax*2];
elseif(method_flag==6)
f_lim_range = [-fmax*2,fmax*2];
end
if plot_flag==1
Power_dB = 10*log10(abs(r1).^2);
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值');
%功率谱估计
yw = fft(r1);
yw=fftshift(yw);
yw = abs(yw).^2/length(yw);
len = length(yw);
f_range = [-fs/2:1/(endT-startT):fs/2];%(0:len-1)/len*fs;
subplot(2,2,2);plot(f_range,10*log10(yw));grid;xlim(f_lim_range);title('周期图法 功率谱 abs()^2/N');xlabel('频率(Hz)');ylabel('功率谱(dB)');
rt = xcorr(r1,r1);
subplot(2,2,3);plot(10*log10(abs(rt)/length(rt)));grid;axis tight;title(['自相关性']);
clear yw;
yw = xcorr(r1,r1);
len = length(yw);
yw = fft(yw);yw=fftshift(yw);
Pyw = abs(yw)/length(yw);%.^2/length(yw);
subplot(2,2,4);
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)');
elseif plot_flag==2
%功率谱估计
fs = 1000/deltaT;
Nfft = 2^4;
while(Nfft)
if (Nfft < length(r1))
Nfft = 2*Nfft;
else
break;
end
end
Nfft = Nfft/2;%让数据长度为2的幂，又不超出采样长度
r2 = r1(1:Nfft);
Power_dB = 20*log10(abs(r2));%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,',method_str]);xlabel('ms');ylabel('dB值');
legend(['E(r^2)=',num2str(10*log10(sum(abs(r2).^2)/length(r2)),'%.2f'),' dB']);
clear Power_dB;
% yw = abs(fftshift(fft(r2))).^2/length(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(f_lim_range);title('周期图法 功率谱');xlabel('频率(Hz)');ylabel('功率谱(dB)');
% yw2=yw;f_less=find(f_range<0);f_more=find(f_range>fmax);yw2([f_less,f_more])=[];
% %size(yw2),fmax*len/fs
% %legend(['(0,fmax)mean=',num2str(10*log10(sum(yw2)*fs/len/fmax),'%.2f'),' dB/Hz']);
% legend(['(0,fmax)mean=',num2str(10*log10(mean(yw2)),'%.2f'),' dB/Hz']);
% clear yw2;
%
% subplot(2,2,3);plot(f_range,10*log10(yw));grid;xlim([-fmax*4,fmax*4]);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(r2))
% 'Welch: 数据总长度应> 1.5*L!'
% return;
% end
%
% while (K>Kmax)
% Lmax = floor(length(r2)*2/L)/2*L; %需要从r2中提取的数据总长度
% K = Lmax*2/L-1; %数据分段数
% if (K <= Kmax)
% if (K==1)
% K = 2;
% L = L/2;
% end
% break;
% else
% L = 2*L;
% end
% end
%
% w_hn = hanning(L);
% Pw = [];
% for k=1:1:K
% Pw(k,:) = (abs(fftshift(fft(w_hn.*r2(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;xlim(f_lim_range);title(['Welch法 功率谱,K=',num2str(K),',L=2\^',num2str(log2(L)),',plot\_flag=',num2str(plot_flag)]);xlabel('频率(Hz)');ylabel('功率谱(dB)');
% f_less=find(f_range<0);f_more=find(f_range>fmax);Pw([f_less,f_more])=[];
% %legend(['(0,fmax)mean=',num2str(10*log10(sum(Pw)*fs/L/fmax),'%.2f'),' dB/Hz']);
% legend(['(0,fmax)mean=',num2str(10*log10(mean(Pw)),'%.2f'),' dB/Hz']);
% clear Pw;
%
% figure;
% subplot(2,2,1);
% [psd_matlab,f_matlab] = periodogram(r2,[],'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));xlim(f_lim_range);grid;
% title(['periodogram(),',num2str(fc), 'MHz,',num2str(v), 'km/h,Max Doppler=',num2str(fmax,'%.2f'),'Hz,']);
% xlabel('Hz');ylabel('dB/Hz');
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;
else
end
y=r1(:);%信道复变量⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -