ray_fading_mea.m

matlab实现的瑞利衰落信道仿真程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 等面积法：本程序只仿真了实部。经典的功率谱密度和自相关函数都和书上相符合。高斯功率谱密度也是书上相符合，但是其自相关函数和书上有出入。
% 还有其自相关函数的均方误差不知道怎么画。还有均值不知道为什么不符合参考模型，难道是因为N不够大？
% N 有限时，符合广义平稳和各态历经，自相关特性与参考模型依然有偏差，因为Ci,n和fi,n是确定的变量。
% N 趋近于无穷时，该模型的统计特性与参考模型相匹配，同时它也满足了广义平稳特性和各态历经性。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;                                  % 清除内存变量和函数
clc;                                    % 清除指令窗
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
startT=0;endT=1; % 单位为s
deltaT=0.001;    % 单位为s

fc=900;         % 载频(MHz)
wc=2*pi*fc;		 % 载波角频率
v=110;            % 绝对时速(km/h)
c=3*10^8;        % 光速(m/s)
fmax = (fc*10^6)*(v*10^3/3600)/c; % Max Doppler Shift (Hz) 最大多普勒频移
wm=2*pi*fmax;							% Maximum shift 最大频移
N = 25; % 模拟均匀到达的rays数 路径数量

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 经典功率谱
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sigma_u = 1;%为rayleigh分布pdf的参数
fd = fmax*sin((1:N)*pi/2/N);  % cos(2*pi*((1:N)/N))*fmax; 
phase = unifrnd(0,2*pi,N,1);  % 相位服从 0-2*pi 均匀分布 ，其中N,1是要产生的随机数的行数和列数
cn = sigma_u*sqrt(2/N);


%计算自相关函数
step = 0;
for t = startT:deltaT:endT/5
	step = step + 1;
	r1(step) = sum(exp(j*(2*pi*fd(:)*t + phase(:))).*cn(:)); % 幅度
    RR(step) = sum((cn(1:end)).^2.*cos(2*pi*fd(1:end)*t)/2);
end
k = (startT:deltaT:endT/5)';
RR1= sigma_u^2*besselj(0,wm*k); % 理论值，当 N 趋近于无穷大时，利用贝塞尔函数画出自相关函数 ，利用P184的（7－1－39）式
r1 = mean(r1(:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 计算功率谱密度
fm4 = [sort(-fd) fd];
L=length(fm4);
Sf = zeros(1,L);

for i=(L/2+1):L
    Sf(i)=cn(i-L/2).^2/4;
    Sf(L-i+1)=Sf(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 计算自相关函数的均方误差
%syms x
%step=0;
%for n=0:5:55
 %   step = step + 1;
 %   taumax = n/2*fmax;
 %   f(step)= 1/taumax*(sigma_u^2*besselj(0,wm*x)-sum((cn(1:end)).^2.*cos(2*pi*fd(1:end)*x)/2))^2;
 %   E(step)=int(f,x,0,taumax);
 %end


figure(1)
plot([startT:deltaT:endT/5],RR,'--r');
title('等距离法 自相关函数,N=25');xlabel('时间(s)');ylabel('r(t)');
hold on;
plot([startT:deltaT:endT/5],RR1,'-b');
legend('--确定型','-理论值');

figure(2)
stem([sort(-fd) fd],Sf);
title('等距离法 功率谱密度,N');xlabel('f/Hz');ylabel('PSD');

figure(3)
plot([startT:deltaT:endT/5],r1,'-r');
title('等距离法 均值,N=25');xlabel('时间(s)');ylabel('均值');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 高斯功率谱
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kc=2*sqrt(2/log(2));
fc1=log(2)^0.5*fmax;
fd2 = kc*fc1*(2*(1:N)-1)/(2*N);    
phase2 = unifrnd(0,2*pi,N,1);   
cn2 = sigma_u*2^0.5*sqrt(erf(kc*log(2)^0.5*(1:N)/N)-erf(kc*log(2)^0.5*((1:N)-1)/N));

%计算自相关函数
step2 = 0;
for t = startT:deltaT*10:endT/5
	step2 = step2 + 1;
    RR2(step2) = sum((cn2(1:end)).^2.*cos(2*pi*fd2(1:end)*t)/2); % 确定型
    RR3(step2) = sigma_u^2*exp(-(pi*fc1*t/log(2)^0.5)^2);   % 理论值
end

% 计算功率谱密度
fm2 = [sort(-fd2) fd2];
L=length(fm2);
Sf2 = zeros(1,L);

for i=(L/2+1):L
    Sf2(i)=cn2(i-L/2).^2/4;
    Sf2(L-i+1)=Sf2(i);
end

figure(4)
plot([startT:deltaT*10:endT/5],RR2,'--r');
title('等距离法 自相关函数,N=25');xlabel('时间(s)');ylabel('r(t)');
hold on;
plot([startT:deltaT*10:endT/5],RR3,'-b');
legend('--确定型','-理论值');

figure(5)
stem([sort(-fd2) fd2],Sf2);
title('等距离法 高斯功率谱密度,N=25');xlabel('f/Hz');ylabel('PSD');