《MIMO-OFDM系统原理、应用及仿真》李莉（实例代码）.docx
《《MIMO-OFDM系统原理、应用及仿真》李莉（实例代码）.docx》由会员分享，可在线阅读，更多相关《《MIMO-OFDM系统原理、应用及仿真》李莉（实例代码）.docx（50页珍藏版）》请在冰点文库上搜索。
2.5仿真实例
实例2-1瑞利分布与莱斯分布
功能：
绘制瑞利分布曲线与莱斯分布曲线
程序名称：
Example2_1.m
程序代码：
clear,clf
N=200000;%产生200000个信道系数供统计使用
level=30;%统计区间被划分的分数。
K_dB=[-40015];%莱斯因子为-40dB、0dB、15dB
gss=['k-*';'k-o';'k-+';'k-^'];%绘制曲线的颜色、线形与标志符号
%瑞利模型
Rayleigh_ch=Ray_model（N）;%调用Ray_model子程序，产生瑞利分布幅度系数
[temp,x]=hist（abs（Rayleigh_ch（1,:
））,level）;%统计数据分布
plot（x,temp,gss（1,:
））
holdon
%莱斯模型
fori=1:
length（K_dB）;%对不同莱斯因子进行信道模型仿真
Rician_ch（i,:
）=Ric_model（K_dB（i）,N）;%调用Ric_model产生莱斯分布幅度系数
[tempx]=hist（abs（Rician_ch（i,:
））,level）;%统计数据分布
plot（x,temp,gss（i+1,:
））
end
xlabel（'x'）,ylabel（'Occurrence'）
legend（'Rayleigh','Rician,K=-40dB','Rician,K=0dB','Rician,K=15dB'）
%瑞利信道模型子程序，子程序程序名称：
Ray_model.m
functionH=Ray_model（L）
%输入参数L:
仿真信道个数，为N=200000
%输出参数H：
返回瑞利信道矩阵
H=（randn（1,L）+j*randn（1,L））/sqrt
（2）;
%产生实部为高斯分布、虚部为高斯分布、包络为瑞利分布的信道系数。
实部功率为1/2，虚部功率为1/2，因
%此该行指令返回单位功率的或称归一化功率的瑞利信道幅度系数。
%莱斯信道模型子程序，子程序程序名称：
Ric_model.m
functionH=Ric_model（K_dB,L）
%输入参数:
K_dB为莱斯因子，L为仿真信道个数
%输出参数H:
返回莱斯信道矩阵
K=10^（K_dB/10）;%将dB值描述的莱斯因子转换为幅度值
H=sqrt（K/（K+1））+sqrt（1/（K+1））*Ray_model（L）;
%产生莱斯信道幅度系数。
莱斯信道模型中包含视距通信，收发之间有直通路径。
程序仿真结果见图2-9。
实例2-2两径信道与指数信道模型
功能：
产生一个两径信道和一个指数衰减的多径信道。
程序名称：
Example2_2
程序代码：
clear,clf
scale=1e-9;%纳秒量级
Ts=10*scale;%抽样时间间隔为10ns，在这个程序中这个量也为指数信道路径间隔
t_rms=30*scale;%RMS时延扩展为30ns
num_ch=10000;%仿真信道个数
%两径信道模型
%产生并绘制了理想的两径信道模型和瑞利分布两径信道模型。
pow_2=[0.50.5];delay_2=[0t_rms*2]/scale;
%给出理想两径信道功率均为0.5，延时为0和60ns
H_2=[Ray_model（num_ch）;Ray_model（num_ch）].'*diag（sqrt（pow_2））;
%产生瑞利两径信道幅度系数。
通过调用子程序Ray_model产生归一化功率的瑞利两径信道幅度系数。
avg_pow_h_2=mean（H_2.*conj（H_2））;
%计算瑞利分布两径信道每一径的平均功率。
在这里可以看到上一条语句中diag（sqrt（pow_2））的作用。
%当通过对幅度系数进行运算计算功率时，sqrt（pow_2）可以使每一径的功率为pow_2，即每一径的功率为0.5。
subplot（121）
stem（delay_2,pow_2,'ko'）,holdon,stem（delay_2,avg_pow_h_2,'k.'）;
xlabel（'Delay[ns]'）,ylabel（'ChannelPower[linear]'）;
title（'2-rayModel'）;
legend（'Ideal','Simulation'）;axis（[-1014000.7]）;
%指数信道模型
%产生并绘制理想的指数信道模型和瑞利分布的指数信道模型。
pow_e=exp_PDP（t_rms,Ts）;%通过调用exp_PDP子程序，计算理想指数信道每一径上的功率。
delay_e=[0:
length（pow_e）-1]*Ts/scale;%计算指数信道每一径的延时，单位为ns
fori=1:
length（pow_e）
H_e（:
i）=Ray_model（num_ch）.'*sqrt（pow_e（i））;
end
%计算瑞利分布的指数信道幅度系数。
通过调用Ray_model产生归一化功率的瑞利分布幅度系数，%sqrt（pow_e（i））的作用类似于diag（sqrt（pow_2））。
avg_pow_h_e=mean（H_e.*conj（H_e））;%计算瑞利分布指数信道的平均功率。
%由于sqrt（pow_e（i））的存在，瑞利分布指数信道每一径的平均功率也为pow_e（i），即与理想指数信道
%每一径功率相同。
subplot（122）
stem（delay_e,pow_e,'ko'）,holdon,stem（delay_e,avg_pow_h_e,'k.'）;
xlabel（'Delay[ns]'）,ylabel（'ChannelPower[linear]'）;
title（'ExponentialModel'）;axis（[-1014000.7]）
legend（'Ideal','Simulation'）
%瑞利信道模型子程序，子程序程序名称：
Ray_model.m
functionH=Ray_model（L）
H=（randn（1,L）+j*randn（1,L））/sqrt
（2）;
%指数信道PDP子程序，子程序名称：
exp_PDP.m
functionPDP=exp_PDP（tau_d,Ts,A_dB,norm_flag）
%输入参数:
%tau_d:
RMS延时扩展，单位为s
%Ts:
抽样时间间隔，在这里也为指数信道路径间隔，单位为s
%A_dB:
最小不可忽略径[dB]
%norm_flag:
标准化标志
%输出参数:
%PDP:
输出指数信道PDP矢量
ifnargin<4,norm_flag=1;end%判断子程序调用参数个数，小于4，则norm_flag=1。
ifnargin<3,A_dB=-20;end%判断子程序调用参数个数，小于4，则A_dB=-20。
%由于主程序中调用该子程序时，只有两个参数，所以上两条语句实际是幅值norm_flag=1和A_dB=-20。
sigma_tau=tau_d;
A=10^（A_dB/10）;
lmax=ceil（-tau_d*log（A）/Ts）;%计算最大路径序号，参见式（2-34）。
%以下参见式（2-36）
ifnorm_flag
p0=（（1-exp（-（lmax+1）*Ts/sigma_tau））/（1-exp（-Ts/sigma_tau）））/30;
elsep0=1/sigma_tau;%计算式（2-37）中的P0
end
%指数信道PDP
l=0:
lmax;
PDP=p0*exp（-l*Ts/sigma_tau）;%参见式（2-37）
程序仿真结果如图2-17所示。
图2-17（a）为理想两径信道和瑞利两径信道的PDP曲线，图2-17（b）为离散指数信道和瑞利指数信道的PDP曲线。
（a）（b）
图2-17两径与指数信道模型
实例2-3IEEE802.11信道PDP与频谱分布
功能：
实现IEEE802.11信道仿真，画出IEEE802.11信道的PDP曲线与频谱图。
程序名称：
Example2_3.m
程序代码：
clear,clf
scale=1e-9;%纳秒量级
Ts=50*scale;%抽样时间间隔，50ns
t_rms=25*scale;%RMS实验扩展，25ns
num_ch=10000;%信道数
N=128;%FFT长度
PDP=IEEE802_11_model（t_rms,Ts）;调用IEEE802_11_model子程序，计算IEEE802.11信道的PDP。
fork=1:
length（PDP）
h（:
k）=Ray_model（num_ch）.*sqrt（PDP（k））;
avg_pow_h（k）=mean（h（:
k）.*conj（h（:
k）））;
end
H=fft（h（1,:
）,N）;
subplot（121）
stem（[0:
length（PDP）-1],PDP,'ko'）,holdon,
stem（[0:
length（PDP）-1],avg_pow_h,'k.'）;
xlabel（'channeltapindex,p'）,ylabel（'AverageChannelPower[linear]'）;
title（'IEEE802.11Model,\sigma_\tau=25ns,T_S=50ns'）;
legend（'Ideal','Simulation'）;axis（[-1701]）;
subplot（122）
plot（[-N/2+1:
N/2]/N/Ts/1e6,10*log10（H.*conj（H））,'k-'）;
xlabel（'Frequency[MHz]'）,ylabel（'Channelpower[dB]'）
title（'Frequencyresponse,\sigma_\tau=25ns,T_S=50ns'）
%IEEE802.11信道模型PDP产生子程序，子程序名称：
IEEE802_11_model.m
functionPDP=IEEE802_11_model（sigma_t,Ts）
%输入参数:
%sigma_t:
RMS延时扩展
%Ts:
抽样时间间隔
%输出参数:
%PDP:
IEEE802.11信道PDP矩阵
lmax=ceil（10*sigma_t/Ts）;%计算最大路径序号，参见式（2-38）
sigma02=（1-exp（-Ts/sigma_t））/（1-exp（-（lmax+1）*Ts/sigma_t））;%参见式（2-41）
l=0:
lmax;PDP=sigma02*exp（-l*Ts/sigma_t）;%参见式（2-40）
仿真结果参见图2-11。
实例2-4滤波白噪声模型
功能：
产生滤波白噪声信道模型，画出信道幅度系数及其包络的概率密度函数与相位的概率密度函数。
程序名称：
Example2_4.m
程序代码：
clear,clf
fm=100;%最大多普勒频率
scale=1e-6;%微秒量级
ts_mu=50;ts=ts_mu*scale;fs=1/ts;%抽样时间与抽样频率
Nd=1e6;%抽样个数
%获得复信道系数
[h,Nfft,Nifft,doppler_coeff]=FWGN_model（fm,fs,Nd）;
subplot（211）
plot（[1:
Nd]*ts,10*log10（abs（h）））%画信道幅度系数
str=sprintf（'Clarke/GanModel,f_m=%d[Hz],T_s=%d[us]',fm,ts_mu）;
title（str）,axis（[00.5-305]）
subplot（223）
hist（abs（h）,50）%画信道包络概率密度函数，包络一维概率密度函数呈瑞利分布。
subplot（224）
hist（angle（h）,50）%画信道相位概率密度函数，相位一维概率密度函数呈均匀分布。
%滤波白噪声信道子程序（Clarke/Gan模型），子程序名称：
FWGN_model.m
function[h,Nfft,Nifft,doppler_coeff]=FWGN_model（fm,fs,N）
%输入参数:
fm为最大多普勒频率，fs为抽样频率，N为抽样个数
%输出参数:
h为复信道系数
Nfft=2^nextpow2（2*fm/fs*N）;
Nifft=ceil（Nfft*fs/（2*fm））;
%产生独立的复高斯随机过程，参见图2-12
GI=randn（1,Nfft）;GQ=randn（1,Nfft）;
%求实信号的FFT，以获得Hermitian对称
CGI=fft（GI）;CGQ=fft（GQ）;%将高斯随机过程转换到频域。
%多普勒谱产生，仿真多普勒滤波器
doppler_coeff=Doppler_spectrum（fm,Nfft）;
%将转换到频域的高斯过程加入到多普勒滤波器，频域内乘积。
f_CGI=CGI.*sqrt（doppler_coeff）;f_CGQ=CGQ.*sqrt（doppler_coeff）;
%补零，使多普勒滤波器输出数据长度为Nifft，以备求ifft。
Filtered_CGI=[f_CGI（1:
Nfft/2）zeros（1,Nifft-Nfft）f_CGI（Nfft/2+1:
Nfft）];
Filtered_CGQ=[f_CGQ（1:
Nfft/2）zeros（1,Nifft-Nfft）f_CGQ（Nfft/2+1:
Nfft）];
%求ifft，将多普勒滤波器输出频域信号转换到时域
hI=ifft（Filtered_CGI）;hQ=ifft（Filtered_CGQ）;
%计算实部的平方加虚部的平方开平方，即呈瑞利分布的包络
rayEnvelope=sqrt（abs（hI）.^2+abs（hQ）.^2）;
%计算包络的均方根值
rayRMS=sqrt（mean（rayEnvelope（1:
N）.*rayEnvelope（1:
N）））;
%图2-12最终输出的Clarke/Gan模型信道冲激响应。
h=complex（real（hI（1:
N））,-real（hQ（1:
N）））/rayRMS;
%多普勒谱子程序，子程序名称：
Doppler_spectrum.m
%对经典多普勒谱的仿真，参见式（2-25）。
functiony=Doppler_spectrum（fd,Nfft）
%输入参数：
fd为最大多普勒频移，Nfft为频域样值点个数。
%输出参数：
y返回多普勒谱
df=2*fd/Nfft;%计算频率间隔
%计算f=0时多普勒谱
f
（1）=0;y
（1）=1.5/（pi*fd）;
%计算其他频率多普勒谱。
计算了多普勒谱从第2个样值到Nfft/2个样值以及第Nfft/2+2个样值到第Nfft个%样值。
fori=2:
Nfft/2,
f（i）=（i-1）*df;%根据频率序号及频率间隔计算频率
y（[iNfft-i+2]）=1.5/（pi*fd*sqrt（1-（f（i）/fd）^2））;
%计算式（2-25），注意randn产生的高斯噪声方差为1。
end
%计算多普勒谱的Nfft/2+1点样值。
用四点样值构成多项式，之后求出多项式在Nfft/2+1点的样值。
这个多%重相当于内插。
nFitPoints=3;kk=[Nfft/2-nFitPoints:
Nfft/2];%四点样值序号为Nfft/2-3到Nfft/2
polyFreq=polyfit（f（kk）,y（kk）,nFitPoints）;%形成多项式
y（（Nfft/2）+1）=polyval（polyFreq,f（Nfft/2）+df）;%从多项式求Nfft/2点值。
仿真结果见图2-13。
3.7仿真实例
实例3-1OFDM信号的产生与解调
功能：
（1）通过对OFDM信号各个子载波赋共轭对称的数据产生一个实OFDM符号；
（2）给OFDM符号加循环前缀与循环后缀；
（3）给OFDM符号加窗。
在程序中加入的是升余弦窗，可以通过改变升余弦窗的滚降系数观察加入不同升余弦窗，对OFDM信号频谱的影响；
（4）信道采用加性高斯白噪声信道。
可以通过改变信噪比改变信道环境，从而在接收端通过误码率或星座图观察信道对OFDM信号传输的影响；
（5）去除循环前缀与循环后缀，对OFDM信号进行解调。
程序名称：
Example3_1.m
程序代码：
clearall;
closeall;
carrier_count=200;
%这个程序中OFDM子载波个数为512，其中400即carrier_count*2为数据符号，其余赋0值。
symbols_per_carrier=20;%每个子载波上的符号数，在这里即为OFDM符号的个数。
bits_per_symbol=4;%OFDM符号的每个子载波上传输的比特数。
4比特通常采用16QAM调制。
IFFT_bin_length=512;%FFT长度，也即一个OFDM符号的子载波的个数。
PrefixRatio=1/4;%循环前缀的比率，即循环前缀与OFDM符号长度的比值，通常在1/6~1/4之间。
GI=PrefixRatio*IFFT_bin_length;%保护间隔的长度，这里为128。
beta=1/32;%升余弦窗的滚降系数。
GIP=beta*（IFFT_bin_length+GI）;%循环后缀的长度，这里为20
SNR=30;%本程序考虑加性高斯白噪声信道，这里信噪比为30dB。
%===============================OFDM信号产生=============================
baseband_out_length=carrier_count*symbols_per_carrier*bits_per_symbol;
%计算传输数据总的比特数，为200*20*4=16000比特。
16000比特的构成为20个OFDM符号，每个OFDM
%符号200个子载波，每个子载波传输4比特信息。
carriers=（1:
carrier_count）+（floor（IFFT_bin_length/4）–floor（carrier_count/2））;
%计算OFDM符号子载波的序号，carriers中存放的序号是29~228。
conjugate_carriers=IFFT_bin_length-carriers+2;
%计算OFDM符号子载波的序号，conjugate_carriers中存放的序号是282~481。
rand（'twister',0）;
baseband_out=round（rand（1,baseband_out_length））;
%产生16000比特待传输的二进制比特流。
这里存放的是发送的二进制信号与后面解调后的二进制信号比
%较，可以计算误码率。
%16QAM调制并绘制星座图
complex_carrier_matrix=qam16（baseband_out）;
%调用子程序qam16进行16QAM调制。
将baseband_out中的二进制比特流，每4比特转换为一个16QAM信
%号，即将二进制比特流每4比特转换为-3-3j、-3+3j、3-3j、3+3j、-1-3j、-1+3j、1-3j、1+3j、%-3-j、-3+j、3-j、3+j、-1-j、-1+j、1-j、1+j中的一个。
转换后complex_carrier_matrix为%1*4000矩阵。
complex_carrier_matrix=reshape…
（complex_carrier_matrix',carrier_count,symbols_per_carrier）';
%转换complex_carrier_matrix中的数据为carrier_count*symbols_per_carrier矩阵，这里为%20*200矩阵。
figure
（1）;
plot（complex_carrier_matrix,'*r'）;%绘制16QAM星座图
axis（[-4,4,-4,4]）;
title（'16QAM调制后星座图'）;
gridon
%IFFT，即进行OFDM调制。
IFFT_modulation=zeros（symbols_per_carrier,IFFT_bin_length）;
%将symbols_per_carrier*IFFT_bin_length矩阵赋0值，这里将20*512矩阵赋0值。
这里512是%IFFT的长度，也是OFDM符号子载波的个数。
IFFT_modulation（:
carriers）=complex_carrier_matrix;
%将20*200的complex_carrier_matrix的数据赋给IFFT_modulation的第29~228列，即给512个子%载波中的29~229个子载波赋值。
IFFT_modulation（:
conjugate_carriers）=conj（complex_carrier_matrix）;
%将20*200的complex_carrier_matrix的数据赋给512个子载波中的第282~481个子载波。
%这段程序构造了512个子载波的OFDM符号，并且各个子载波上的数据是共轭对称的。
这样做的目的是经过%IFFT后形成的OFDM符号均为实数。
另外，在512个子载波中，仅有400个子载波为数据，其余为0值。
相%当于补零，补零的目的是通常IFFT的长度应该为2的整数次幂。
signal_after_IFFT=ifft（IFFT_modulation,IFFT_bin_length,2）;%IFFT实现OFDM调制。
time_wave_matrix=signal_after_IFFT;%
figure
（2）;
plot（0:
IFFT_bin_length-1,time_wave_matrix（2,:
））;%画一个OFDM信号的时域表现
axis（[0,512,-0.4,0.4]）;
gridon;
ylabel（'Amplitude'）;
xlabel（'Time'）;
title（'OF
升级会员 