引言
最近完成老师给的作业,题目如下:
无线信道中的多普勒谱有一种经典谱(classic spectrum)。请用Simulink或者m语言,产生一条单径瑞利信道,其多普勒谱为经典谱,其中移动速率为120km/h。
在查阅了多普勒经典谱的相关知识后,首先接触到的就是Clarke模型、Jakes模型这两个模型。在此先阐述他们的区别:
Clarke模型是一个多径衰落的数学模型
Jakes模型是是一个仿真模型,(也就是说是clarke模型的一个具体实现形式)它产生的信号是广义平稳的,并且能够较好的吻合Clarke模型中的统计特性
经典谱的产生(Clarke模型,即理论推导)
所谓的经典谱,最早是由Clarke于1968年提出来的。这个模型假设:
存在N个有任意相位的平面波,每个平面波以任意的方向达到接收机,且平均功率相同
也就是说,平面波达到接收机的仰角为零,而水平方位角(暂记为θ)在0到2 π \pi π弧度之间均匀连续分布。
我们任选其中一束平面波来分析:
假设,发射信号为:
x
(
t
)
x(t)
x(t)
则通频带发射信号:
x
~
(
t
)
=
R
e
[
x
(
t
)
e
j
2
π
f
c
t
]
\tilde{x}(t) = {\mathop{\rm Re}\nolimits} \left[ {x(t){e^{j2\pi {f_c}t}}} \right]
x~(t)=Re[x(t)ej2πfct]
由于移动台的运动,到达接收机的所有平面波都会经历多普勒频移。则通过
I
I
I条传播路径的散射信道后(不同路径多普勒频移不同),通频带接收信号可表示为:
y
~
(
t
)
=
R
e
[
∑
i
=
1
I
C
i
e
j
2
π
(
f
c
+
f
i
)
(
t
−
τ
i
)
x
(
t
−
τ
i
)
]
=
R
e
[
y
(
t
)
e
j
2
π
f
c
t
]
\tilde{y}(t) = {\mathop{\rm Re}\nolimits} \left[ {\sum\limits_{i = 1}^I {{C_i}} {e^{j2\pi ({f_c} + {f_i})(t - {\tau _i})}}x(t - {\tau _i})} \right]= {\mathop{\rm Re}\nolimits} \left[ {y(t){e^{j2\pi {f_c}t}}} \right]
y~(t)=Re[i=1∑ICiej2π(fc+fi)(t−τi)x(t−τi)]=Re[y(t)ej2πfct]
(注:
C
i
C_i
Ci,
τ
i
\tau_i
τi,
f
i
f_i
fi分别表示第i条路径的信道增益,时延和多普勒频移。其中
f
i
=
f
m
c
o
s
θ
=
v
λ
c
o
s
θ
f_i=f_mcos\theta=\frac{v}{\lambda}cos\theta
fi=fmcosθ=λvcosθ,
f
m
f_m
fm为最大多普勒频移)
基带接收信号为:
y
(
t
)
=
∑
i
=
1
I
C
i
e
−
j
ϕ
i
(
t
)
x
(
t
−
τ
i
)
y(t) = \sum\limits_{i = 1}^I {{C_i}{e^{ - j{\phi _i}(t)}}} x(t - {\tau _i})
y(t)=i=1∑ICie−jϕi(t)x(t−τi)
其中
ϕ
i
(
t
)
=
2
π
{
(
f
c
+
f
i
)
τ
i
−
f
i
t
i
}
{\phi _i(t)}=2\pi\{(f_c+f_i)\tau_i-f_it_i\}
ϕi(t)=2π{(fc+fi)τi−fiti}
所以尽管发射频率为
f
c
f_c
fc,但是接收到的信号频谱却扩展到了
f
c
−
f
m
f_c-f_m
fc−fm到
f
c
+
f
m
f_c+f_m
fc+fm范围内,这就是多普勒频展。
根据输入输出,可以把信道建模为一个线性时变滤波器,脉冲响应如下:
h
(
t
,
τ
)
=
∑
i
=
1
I
C
i
e
−
j
ϕ
i
(
t
)
δ
(
t
−
τ
i
)
≈
h
(
t
)
δ
(
t
−
τ
^
)
h(t,\tau ) = \sum\limits_{i = 1}^I {{C_i}{e^{ - j{\phi _i}(t)}}} \delta (t - {\tau _i}) \approx h(t)\delta (t - \hat{\tau} )
h(t,τ)=i=1∑ICie−jϕi(t)δ(t−τi)≈h(t)δ(t−τ^)
其中
h
(
t
)
=
∑
i
=
1
I
C
i
e
−
j
ϕ
i
(
t
)
h(t)=\sum\limits_{i = 1}^I {{C_i}{e^{ - j{\phi _i}(t)}}}
h(t)=i=1∑ICie−jϕi(t)。
将此表达式回带进通频带接收信号,假设
x
(
t
)
=
1
x(t)=1
x(t)=1,可得通频带接收信号表达式为:
y
~
(
t
)
=
R
e
[
y
(
t
)
e
j
2
π
f
c
t
]
=
R
e
[
{
h
I
(
t
)
+
j
h
Q
(
t
)
}
e
j
2
π
f
c
t
]
=
h
I
(
t
)
cos
2
π
f
c
t
−
h
Q
(
t
)
sin
2
π
f
c
t
\begin{array}{c} \tilde{y}(t) = {\mathop{\rm Re}\nolimits} \left[ {y(t){e^{j2\pi {f_c}t}}} \right]\\ = {\mathop{\rm Re}\nolimits} \left[ {\{ {h_I}(t) + j{h_Q}(t)\} {e^{j2\pi {f_c}t}}} \right]\\ = {h_I}(t)\cos 2\pi {f_c}t - {h_Q}(t)\sin 2\pi {f_c}t \end{array}
y~(t)=Re[y(t)ej2πfct]=Re[{hI(t)+jhQ(t)}ej2πfct]=hI(t)cos2πfct−hQ(t)sin2πfct
其中
h
I
(
t
)
h_I(t)
hI(t)和
h
Q
(
t
)
h_Q(t)
hQ(t)是
h
(
t
)
h(t)
h(t)的同相和正交分量。
h
I
(
t
)
=
∑
i
=
1
I
C
i
cos
ϕ
i
(
t
)
h_I(t)=\sum\limits_{i = 1}^I {{C_i}\cos {\phi _i}(t)}
hI(t)=i=1∑ICicosϕi(t)
h
Q
(
t
)
=
∑
i
=
1
I
C
i
sin
ϕ
i
(
t
)
h_Q(t)=\sum\limits_{i = 1}^I {{C_i}\sin {\phi _i}(t)}
hQ(t)=i=1∑ICisinϕi(t)
只要
I
I
I足够大,根据中心极限定理我们可知,
h
I
(
t
)
h_I(t)
hI(t)和
h
Q
(
t
)
h_Q(t)
hQ(t)会趋近高斯随机变量,此时接收信号的幅度
∣
y
~
(
t
)
∣
=
h
I
2
(
t
)
+
h
Q
2
(
t
)
|\tilde{y}(t)| = \sqrt {{h_I}^2(t) + {h_Q}^2(t)}
∣y~(t)∣=hI2(t)+hQ2(t)服从瑞利分布
对
y
~
(
t
)
\tilde{y}(t)
y~(t)的自相关函数做傅里叶变换,得到衰落过程的功率谱密度
这就是经典多普勒谱。话不多说,让我们亲自见识一下
经典谱的产生(Jakes模型,即仿真实践)
单径瑞利信道,是信道的一种窄带衰落模型,在这种衰落模型中,信道时延扩展相对很小,多径分量不可分辨,信道为平坦衰落信道。可以用小信号模型来进行建模和拟合。
根据理论模型我们可知:复高斯信号的模服从瑞利分布,因此就需要产生一个复高斯信号,虚部实部都服从高斯分布。
那么我们如何生成这样的实部虚部呢?我们采用正弦波叠加法,理论依据是中心极限定理:无穷多个独立同分布信号的叠加服从高斯分布。
假设我们有N个平面波,原始信号的频率是
ω
d
=
2
π
f
m
=
2
π
v
λ
\omega_d=2\pi f_m=2\pi \frac{v}{\lambda}
ωd=2πfm=2πλv,相位是
ϕ
N
=
0
\phi_N=0
ϕN=0
定义 N 0 = ( N / 2 − 1 ) / 2 N_0=(N/2-1)/2 N0=(N/2−1)/2,其中N/2被限定为一个奇数。
如上图所示,将
N
0
N_0
N0个频率为
ω
n
\omega_n
ωn,相位为
ϕ
n
\phi_n
ϕn(均匀分布)的复振荡器的输出求和,然后和频率为
ω
d
=
2
π
f
m
\omega_d=2\pi f_m
ωd=2πfm的复振荡器的输出相加。
(每个复振荡器的频率为
ω
n
=
ω
d
c
o
s
θ
n
\omega_n=\omega_dcos\theta_n
ωn=ωdcosθn,
θ
n
=
2
π
n
N
\theta_n=\frac{2\pi n}{N}
θn=N2πn表示第n个平面波的到达角度。
每个复振荡器的相位为
ϕ
n
=
π
n
N
0
+
1
\phi_n=\frac{\pi n}{N_0+1}
ϕn=N0+1πn。
由此求和得到实部
h
I
(
t
)
h_I(t)
hI(t)和虚部
h
Q
(
t
)
h_Q(t)
hQ(t)(为便于编程,下面计算都穿插了矩阵运算)
h
I
(
t
)
=
2
∑
n
=
1
N
0
(
cos
ϕ
n
cos
ω
n
t
)
+
2
cos
ϕ
N
cos
ω
d
t
=
[
cos
ϕ
1
cos
ϕ
2
⋯
cos
ϕ
N
0
cos
ϕ
N
]
⋅
[
2
cos
ω
1
t
2
cos
ω
2
t
⋮
2
cos
ω
N
0
t
2
cos
ω
d
t
]
h{}_I(t) = 2\sum\limits_{n = 1}^{{N_0}} {\left( {\cos {\phi _n}\cos {\omega _n}t} \right) + \sqrt 2 \cos } {\phi _N}\cos {\omega _d}t\\= \begin{bmatrix} {\cos {\phi _1}}&{\cos {\phi _2}}& \cdots &{\cos {\phi _{{N_0}}}}&{\cos {\phi _N}} \end{bmatrix}\cdot \begin{bmatrix} {2\cos {\omega _1}t}\\ {2\cos {\omega _2}t}\\ \vdots \\ {2\cos {\omega _{{N_0}}}t}\\ {2\cos {\omega _d}t} \end{bmatrix}
hI(t)=2n=1∑N0(cosϕncosωnt)+2cosϕNcosωdt=[cosϕ1cosϕ2⋯cosϕN0cosϕN]⋅⎣⎢⎢⎢⎢⎢⎡2cosω1t2cosω2t⋮2cosωN0t2cosωdt⎦⎥⎥⎥⎥⎥⎤
h Q ( t ) = 2 ∑ n = 1 N 0 ( sin ϕ n cos ω n t ) + 2 sin ϕ N cos ω d t = [ sin ϕ 1 sin ϕ 2 ⋯ sin ϕ N 0 sin ϕ N ] ⋅ [ 2 cos ω 1 t 2 cos ω 2 t ⋮ 2 cos ω N 0 t 2 cos ω d t ] h{}_Q(t) = 2\sum\limits_{n = 1}^{{N_0}} {\left( {\sin {\phi _n}\cos {\omega _n}t} \right) + \sqrt 2 \sin } {\phi _N}\cos {\omega _d}t\\= \begin{bmatrix} {\sin {\phi _1}}&{\sin {\phi _2}}& \cdots &{\sin {\phi _{{N_0}}}}&{\sin {\phi _N}} \end{bmatrix}\cdot \begin{bmatrix} {2\cos {\omega _1}t}\\ {2\cos {\omega _2}t}\\ \vdots \\ {2\cos {\omega _{{N_0}}}t}\\ {2\cos {\omega _d}t} \end{bmatrix} hQ(t)=2n=1∑N0(sinϕncosωnt)+2sinϕNcosωdt=[sinϕ1sinϕ2⋯sinϕN0sinϕN]⋅⎣⎢⎢⎢⎢⎢⎡2cosω1t2cosω2t⋮2cosωN0t2cosωdt⎦⎥⎥⎥⎥⎥⎤
Jakes模型的复输出可以表示为
h
(
t
)
=
E
0
2
N
0
+
1
{
h
I
(
t
)
+
j
h
Q
(
t
)
}
h(t) = \frac{{{E_0}}}{{\sqrt {2{N_0} + 1} }}\{ {h_I}(t) + j{h_Q}(t)\}
h(t)=2N0+1E0{hI(t)+jhQ(t)}
其中
E
0
E_0
E0为衰落信道的平均幅度。
h
I
(
t
)
+
j
h
Q
(
t
)
{h_I}(t) + j{h_Q}(t)
hI(t)+jhQ(t)带入矩阵运算可表示为
h
I
(
t
)
+
j
h
Q
(
t
)
=
[
cos
ϕ
1
+
j
sin
ϕ
1
cos
ϕ
2
+
j
sin
ϕ
2
⋯
cos
ϕ
N
0
+
j
sin
ϕ
N
0
cos
ϕ
N
+
j
sin
ϕ
N
]
⋅
[
2
cos
ω
1
t
2
cos
ω
2
t
⋮
2
cos
ω
N
0
t
2
cos
ω
d
t
]
=
[
e
j
ϕ
1
e
j
ϕ
1
⋯
e
j
ϕ
N
0
e
j
ϕ
N
]
⋅
[
2
cos
ω
1
t
2
cos
ω
2
t
⋮
2
cos
ω
N
0
t
2
cos
ω
d
t
]
=
2
∑
n
=
1
N
0
e
j
ϕ
n
cos
ω
n
t
+
2
e
j
ϕ
N
cos
ω
d
t
{h_I}(t) + j{h_Q}(t)= \begin{bmatrix} {\cos {\phi _1}+j\sin {\phi _1}}&{\cos {\phi _2}+j\sin {\phi _2}}& \cdots &{\cos {\phi _{N_0}}+j\sin {\phi _{{N_0}}}}&{\cos {\phi _{N}}+j\sin {\phi _N}} \end{bmatrix}\cdot \begin{bmatrix} {2\cos {\omega _1}t}\\ {2\cos {\omega _2}t}\\ \vdots \\ {2\cos {\omega _{{N_0}}}t}\\ {2\cos {\omega _d}t} \end{bmatrix}\\= \begin{bmatrix} {e^{j \phi_1}}&{e^{j \phi_1}}& \cdots &{e^{j \phi_{N_0}}}&{e^{j \phi_N}} \end{bmatrix}\cdot \begin{bmatrix} {2\cos {\omega _1}t}\\ {2\cos {\omega _2}t}\\ \vdots \\ {2\cos {\omega _{{N_0}}}t}\\ {2\cos {\omega _d}t} \end{bmatrix}\\=2\sum\limits_{n = 1}^{{N_0}} {e^{j\phi_n}\cos {\omega _n}t} + \sqrt 2 e^{j\phi_N}\cos {\omega _d}t
hI(t)+jhQ(t)=[cosϕ1+jsinϕ1cosϕ2+jsinϕ2⋯cosϕN0+jsinϕN0cosϕN+jsinϕN]⋅⎣⎢⎢⎢⎢⎢⎡2cosω1t2cosω2t⋮2cosωN0t2cosωdt⎦⎥⎥⎥⎥⎥⎤=[ejϕ1ejϕ1⋯ejϕN0ejϕN]⋅⎣⎢⎢⎢⎢⎢⎡2cosω1t2cosω2t⋮2cosωN0t2cosωdt⎦⎥⎥⎥⎥⎥⎤=2n=1∑N0ejϕncosωnt+2ejϕNcosωdt
此外,经过多普勒频移的正弦数
N
0
N_0
N0必须足够大,以便衰落信道的振幅能够近似服从瑞利分布。一般取
N
0
N_0
N0=8就足够大了。
得到结果后,我们用图表的形式展现出来(时域特性、幅度谱、相位谱、自相关函数、功率谱):
代码如下,有需要自取,烦请标明引用出处.
记得帮忙一键三连,或者点个赞也行
主函数
close all, clear all,clc
%参数
fd= 926; %多普勒频率
Ts= 1e-6;%采样周期
M= 2^12;
t=[0:M-1]*Ts;
f=[-M/2:M/2-1]/(M*Ts*fd);
Ns = 50000;
t_state=0; %%最终时间
%信道生成
[h,t_state]= Jakes_Flat(fd,Ts,Ns,t_state,1,0);
subplot(311)
plot([1:Ns]*Ts,10*log10(abs(h)))
title(['Jakes Model, f_d=',num2str(fd),'Hz,T_s=',num2str(Ts),'s']);
axis([0 Ns*Ts -20 10])
xlabel('time[s]'),ylabel('Magnitude[dB]');
subplot(323)
hist(abs(h), 50)
title(['Jakes Model, f_d=',num2str(fd),'Hz,T_s=',num2str(Ts),'s']);xlabel('Magnitude'), ylabel('Occasions')
subplot(324)
hist(angle(h), 50)
title(['Jakes Model,f_d=',num2str(fd),'Hz,T_s=',num2str(Ts),'s']);xlabel('Phase[rad])'),ylabel('Occasions')
%信道的自相关函数
temp = zeros(2,Ns);
for i =1: Ns
j= i:Ns;
temp(1:2,j-i+1)= temp(1:2,j-i+1)+[h(i)'*h(j);ones(1,Ns-i+1)];
end
for k=1:M
Simulatcd_corr(k) = real(temp(1,k))./temp(2,k);
end
Classical_corr = besselj(0,2*pi*fd*t);
%自相关函数的傅里叶变换
Classical_Y= fftshift(fft(Classical_corr));
Simulated_Y= fftshift(fft(Simulatcd_corr));
subplot(325)
plot(t,abs(Classical_corr),'k-',t,abs(Simulatcd_corr),'r:')
title(['Autocorrelation, f d=',num2str(fd),'Hz'])
grid on
xlabel('delay \tau[s]'), ylabel('Correlation')
legend('Classical','Simulated')
subplot(326)
plot(f,abs(Classical_Y),'k-', f,abs(Simulated_Y),'r:')
title(['Doppler Spectrum,f_d=',num2str(fd),'Hz'])
axis([-1 1 0 600]),
xlabel('f/f_d'), ylabel('Magnitude')
legend('Classical','Simulated')
信道模型函数
function [h,tf]= Jakes_Flat(fd,Ts,Ns,t0,EO,phi_N) %Jakes 信道模型
%输入:
% fd:多普勒频率
% Ts:采样周期
% Ns:采样点数
% t0:初始时间
% EO:信道功率
% phi_N:具有最大多普勒频率正弦信号的初始相位
%输出:
% h:复衰落向量
% t_state:当前时刻
if nargin<6,phi_N=0;end
if nargin<5,EO=1; end
if nargin<4,t0=0;end
if nargin<3
error('More arguments are needed for Jakes_Flat()');
end
NO = 8;
N=4*NO+2;
wd=2*pi*fd;%最大多普勒频率[rad]
t=t0+[0:Ns-1]*Ts;%时间向量
tf = t(end)+Ts;%最终时间
coswt =[sqrt(2)*cos(0)*cos(wd*t); 2*cos(wd*cos(2*pi*([1:NO]')/N)*t)];
h= EO/sqrt(2*NO+1)*exp(j*[phi_N pi*([1:NO])/(NO+1)])*coswt;
end
本文及程序参考自:《MIMO-OFDM Wireless Communications with MATLAB》