基于Matlab----多径衰落信道

一、题目

多径信道仿真

二、仿真要求

通过一个简单的模拟程序来说明多径衰落信道的特点，针对影响信道的两个重要参数2径，移动台速度来说明相干带宽、相干时间的定义。

三、主要代码

本次实验有三种模型可供选择：Clarke模型、Jakes模型和Zheng模型。我选用的是Jakes模型，但我也会提供另外两种模型的代码。
程序代码：多径衰落信道matlab仿真代码

3.1、Jakes模型

Jakes模型代码：

clear;clc
fd=926;     % 最大Doppler
ts_mu=50;   % 持续时间
scale=1e-6;
ts=ts_mu*scale; %采样时间
fs=1/ts;    % sample frequency
Ns=1e5;     % 采样数

M=2^12;
t=(0:Ns-1)*ts;
tt=(0:M-1)*ts;
ff=[-M/2:M/2-1]/(M*ts*fd);
temp=zeros(3,Ns);

%generat channel information
for ii=1:50
    t_state=0;  % 开始采样时刻
    [h,t_state]=Jakes_model(fd,ts,Ns,t_state,1,0);  % generate channel
    yy=xcorr(h);
    yy_cs=xcorr(real(h),imag(h));
    temp(1,:)=temp(1,:)+yy(Ns:length(yy));
    temp(3,:)=temp(3,:)+yy_cs(Ns:length(yy_cs));
end

figure(1);
subplot(211)
plot((1:Ns)*ts,10*log10(abs(h)));
%axis([0 0.5 -20 5]);
axis([0 0.01 -20 10]);
xlabel('时间/s');ylabel('幅度/db');title('信道时域特性');
str=sprintf('channel model by Jakes with fm=%d[Hz],Ts=%d[us]',fd,ts_mu);
title(str);
% 信道包络
subplot(223)
[f,xi]=ksdensity(abs(h));
plot(xi,f);
cs2=var(h)/2;  %公式乘上了系数sqrt(2)
r=linspace(0,2,1000);
fx2=r./(cs2).*exp(-r.^2/2/(cs2));
hold on;plot(r,fx2,'r:');hold off;
xlabel('幅度');ylabel('统计次数');title('幅度分布');axis([0 2.5 0 1.0]);
% 信道相位
subplot(224)
[f,xi]=ksdensity(angle(h));
plot(xi,f);
hold on;plot([-pi pi],[1/2/pi 1/2/pi],'r:');hold off;
xlabel('相位/rad');ylabel('统计次数');title('相位分布');axis([-pi pi 0 0.2]);


%%自相关函数和功率谱密度
temp(1,:)=temp(1,:)/50;
temp(3,:)=temp(3,:)/50;
%用于归一化
yyy=xcorr(ones(1,Ns));
temp(2,:)=yyy(Ns:length(yy));
%自相关和分量互相关
for k=1:M
    simulated_corr(k)=real(temp(1,k))/temp(2,k);
    simulated_corr_cs(k)=real(temp(3,k))/temp(2,k);
end
classical_corr=besselj(0,2*pi*fd*tt);

%%画图：自相关函数和互相关函数
figure(2);
subplot(211)
plot(tt,classical_corr,'k-',tt,simulated_corr,'r-');
title('自相关函数'); grid on; 
xlabel('时间差/s');ylabel('相关系数');axis([0 0.004 -0.5 1]);
legend('理想特性','仿真结果');
subplot(212)
plot(tt,simulated_corr_cs,'r-');
title('互相关函数'); grid on; 
xlabel('时间差/s');ylabel('相关系数');axis([0 0.004 -0.5 0.5]);

3.2、Clarke模型

Clarke模型代码：

clear;clc
fd=926;     % 最大Doppler
ts=5e-5;    % 仿真最小时间
fs=1/ts;    % sample frequency
Ns=1e5;     % 仿真次数
M=2^12;     % 做FFT的总点数，可以使用自相关函数通过FFT计算功率谱密度，请自行完成
t=(0:Ns-1)*ts;% 仿真总时间，作图使用
tt=(0:M-1)*ts;% 仿真总时间，作图使用
ff=(-M/2:M/2-1)/(M*ts*fd);% 频率归一化，功率谱密度作图使用
temp=zeros(3,Ns); % 3*Ns的矩阵
                  % temp（1，：）存储仿真结果的自相关函数
                  % temp（2，：）存储ones（1，Ns）的自相关函数，用于归一化
                  % temp（3，：）存储仿真结果其分量的互相关函数
                  
%generat channel information
[h,Nfft,Nifft,Doppler_coeff]=Clarke_model(fd,fs,Ns); % 归一化的Clarke模型，功率为1

%自相关函数，分量的互相关函数，为了结果更准确，这里50次仿真取平均值
for ii=1:50
    [h,Nfft,Nifft,Doppler_coeff]=Clarke_model(fd,fs,Ns);
    yy=xcorr(h);
    yy_cs=xcorr(real(h),imag(h));
    temp(1,:)=temp(1,:)+yy(Ns:length(yy));
    temp(3,:)=temp(3,:)+yy_cs(Ns:length(yy_cs));
end

%画图
%Clarke模型复包络的幅度增益、幅度分布、相位分布
figure(1);
subplot(211)
plot((1:Ns)*ts,10*log10(abs(h)));
axis([0 0.01 -20 10]);
xlabel('时间/s');ylabel('幅度/db');title('信道时域特性');
str=sprintf('channel model by Clarke with fm=%dHz,Ts=%dus',fd,ts*1e6);
title(str);
% 幅度分布
subplot(223)
[f,xi]=ksdensity(abs(h));
plot(xi,f);
cs2=var(h)/2;  % 瑞利分布的方差是实部或虚部的方差，是h方差的一半
r=linspace(0,2,1000);
fx2=r./(cs2).*exp(-r.^2/2/(cs2));
hold on;plot(r,fx2,'r:');hold off;
xlabel('幅度');ylabel('统计次数');title('幅度分布');axis([0 2.5 0 1.0]);
% 相位分布
subplot(224)
[f,xi]=ksdensity(angle(h));
plot(xi,f);
hold on;plot([-pi pi],[1/2/pi 1/2/pi],'r:');hold off;
xlabel('相位/rad');ylabel('统计次数');title('相位分布');axis([-pi pi 0 0.2]);


%%自相关函数和互相关函数计算，功率谱密度请自行编程
temp(1,:)=temp(1,:)/50;
temp(3,:)=temp(3,:)/50;
%用于归一化
yyy=xcorr(ones(1,Ns));
temp(2,:)=yyy(Ns:length(yy));
%自相关函数和分量互相关函数，取M个点，进行归一化
for k=1:M
    simulated_corr(k)=real(temp(1,k))/temp(2,k);
    simulated_corr_cs(k)=real(temp(3,k))/temp(2,k);
end
classical_corr=besselj(0,2*pi*fd*tt);

%%画图：自相关函数和互相关函数
figure(2);
subplot(211)
plot(tt,simulated_corr,'k-');
title('自相关函数'); grid on; 
xlabel('时间差/s');ylabel('相关系数');axis([0 0.004 -0.5 1]);
legend('仿真结果');
subplot(212)
plot(tt,simulated_corr_cs,'k-');
title('互相关函数'); grid on; 
xlabel('时间差/s');ylabel('相关系数');axis([0 0.004 -0.2 0.2]);

3.3、Zheng模型

Zheng模型代码：

clear;clc;
fd=926;
ts=5e-5;
Ns=1e5;
M=2^12;

t=(0:Ns-1)*ts;
tt=(0:M-1)*ts;
ff=[-M/2:M/2-1]/(M*ts*fd);
temp=zeros(3,Ns);

for ii=1:50
    h=z_rayleigh(8,fd,t);
    yy=xcorr(h);
    yy_cs=xcorr(real(h),imag(h));
    temp(1,:)=temp(1,:)+yy(Ns:length(yy));
    temp(3,:)=temp(3,:)+yy_cs(Ns:length(yy_cs));
end

%%自相关函数和功率谱密度
temp(1,:)=temp(1,:)/50;
temp(3,:)=temp(3,:)/50;
%用于归一化
yyy=xcorr(ones(1,Ns));
temp(2,:)=yyy(Ns:length(yy));
%自相关和分量互相关
for k=1:M
    simulated_corr(k)=real(temp(1,k))/temp(2,k);
    simulated_corr_cs(k)=real(temp(3,k))/temp(2,k);
end
classical_corr=besselj(0,2*pi*fd*tt);
%功率谱密度，不再作图
classical_y=fftshift(fft(classical_corr));
simulated_y=fftshift(fft(simulated_corr));

%%画图：自相关函数和互相关函数
figure(2);
subplot(211)
plot(tt,classical_corr,'k-',tt,simulated_corr,'r-');
title('自相关函数'); grid on; 
xlabel('时间差/s');ylabel('相关系数');axis([0 0.004 -0.5 1]);
legend('理想特性','仿真结果');
subplot(212)
plot(tt,simulated_corr_cs,'r-');
title('互相关函数'); grid on; 
xlabel('时间差/s');ylabel('相关系数');axis([0 0.004 -0.2 0.2]);



%幅度和相位分布函数
[f,xi]=ksdensity(abs(h));
[f2,xii]=ksdensity(angle(h));
%rayleigh分布的pdf
cs2=var(h)/2;  %方差为实部或虚部的方差
r=linspace(0,2,1000);
fx2=r./(cs2).*exp(-r.^2/2/(cs2));

%%画图：幅度、幅度分布、相位分布
figure(1);
%幅度增益
subplot(211);
plot(t,10*log10(abs(h)));
ylabel('幅度/dB');xlabel('时间/s');axis([0 0.01 -20 10]);
title('改进Jakes模型，fd=926Hz，ts=50us');
subplot(223);
plot(xi,f);
hold on;plot(r,fx2,'r:');hold off;
ylabel('幅度分布');xlabel('幅度');title('幅度分布');axis([0 2.5 0 1]);
subplot(224);
plot(xii,f2);
hold on;plot([-pi pi],[1/2/pi 1/2/pi],'r:');hold off;
ylabel('相位分布');xlabel('相位/rad');title('相位分布');axis([-pi pi 0 0.2]);

四、仿真结果：

Jakes模型仿真结果：