引言

最近“无线通信”课程的仿真实验，要求我们进行针对几种常见的无线通信信道进行仿真实现。接下来，本文将从瑞利衰落信道的理论、相关模型出发，分析基于瑞利衰落信道的多径衰落信道。

由于现在是大三，将来有志于升学，所以这些通信方面的知识算是地基了。在这里算是一个实验心得笔记叭~

以下是我的一些分享。如有不正确的地方，还请指正。

小尺度衰落理论

在梳理瑞利信道的知识点之前，首先需要回顾一下基础中的基础——小尺度衰落。

在“无线通信”这门课中，比较难的就是这块小尺度衰落。刚刚开课那几天，老师每节课必提一嘴“小尺度衰落这一块是比较难的”，“历届讲到二维傅里叶变换的时候，一遍过去没人能听懂的！”——闲话少叙，那就进入理论梳理吧~

无线信道的小尺度衰弱特征可以分为三大类：

• 一类是由于多径传播导致短时间内幅度衰落；
• 一类是由于多径的时延扩展引起时间色散导致的信道衰弱；
• 一类是由于多普勒扩展引起频率色散导致的信道衰弱。

当发射信号通过无线信道传播时，信号参数和信道时间色散与频率色散参数之间的关系决定了发射信号所经历的小尺度衰弱类型。

多径时延扩展产生的衰落效应

多径特性引起的时间色散，导致发送的信号产生平坦衰落或频率选择性衰落。

平坦衰落

如果移动无线信道带宽大于发送信号的带宽，且在带宽范围内有恒定增益及线性相位，则接收信号就会经历平坦衰落过程。这种衰落是最常见的一种。在平坦衰落情况下，信道的多径结构使发送信号的频谐特性在接收机内仍能保持不变。然而，由于多径导致信道增益的起伏，使接收信号的强度会随着时间变化。平坦衰落信道的特性如图所示。

在平坦衰落信道中，发送信号带宽的倒数远大于信道的多径时延扩展，$h_b(t，\tau )$可近似认为无附加时延（即τ=0的单一δ函数）。平坦衰落信道即幅度变化信道，有时看成窄带信道，这是由于信号带宽比平坦衰落信道带宽窄得多。典型的平坦衰落信道会引起深度衰落，因此在深度衰落期间需要增加20dB或30dB的发送功率，以获得较低的比特误码率。

其中，确切来说，如果相干带宽定义为频率相关函数大于0.9 的某特定带宽，则相干带宽近似为
$$B_c=1/50{\sigma }_{\tau }$$
如果将定义放宽至相关函数值大于0.5 ，则相干带宽近似为：
$$B_c=1/5{\sigma }_{\tau }$$
经历平坦衰落的条件为$B_s\ll B_{\tau }$或$T_s\gg {\sigma }_{\tau }$，近似估计时取$T_s\ge 10{\sigma }_{\tau }$。

频率选择性衰落

如果信道具有恒定增益和线性相位的带宽范围小于发送信号带宽，则该信道特性会导致接收信号产生选择性衰落。在这种情况下，信道冲激响应具有多径时延扩展，其值大于发送信号波形带宽的倒数。此时，接收信号中包含经历了衰减和时延的发送信号波形的多径波， 多径信道引入的符号间干扰(ISI)明显，因而，产生接收信号失真。频率选择性衰落信道的特性如图所示。

频率选择性衰落信道的建模比平坦衰落信道的建模更困难，因为必须对每一个多径信号建模，且必须把信道视作一个线性滤波器。为此要进行宽带多径测量，并在此基础上进行建模。

通常频率选择性衰弱信道的判定条件为$B_s>B_{\tau }$或$T_s<{\sigma }_{\tau }$，近似估计时取$T_s<10{\sigma }_{\tau }$。

多普勒扩展引起的衰落效应

快衰落

根据发送信号与信道变化快慢程度的比较，信道可分为快衰落信道和慢衰落信道。在快衰落信道中，信道冲激响应在符号周期内变化很快。即信道的相干时间比发送信号的信号周期短。由于多普勒扩展引起频率色散(也称为时间选择性衰落)，从而导致信号失真。从频域可看出，信号失真随发送信号带宽的多普勒扩展的增加而加剧。因此，信号经历快衰落的条件是：
$$B_d>W_s$$$$(\Delta t{)}_c<T_s$$
其中$B_d$为多普勒扩展，$(\Delta t{)}_c$为信道相干时间，$(\Delta t{)}_c\approx 1/B_d$ ，$W_s$为信号带宽，$T_s$为信号周期，$T_s=1/W_s$ 。

慢衰落

在慢衰落信道中，信道冲激响应变化率比发送的基带信号S(t)变化率低得多，因此可假设在一个或若干个带宽倒数间隔内，信道均为静态信道。在频域中，这意味着信道的多普勒扩展比基带信号带宽小的多。所以信号经历慢衰落的条件是：
$$B_d\ll W_s$$$$(\Delta t{)}_c\gg T_s$$
显然，移动台的速度(或信道路径中物体的速度)及基带信号发送速率，决定了信号是经历快衰落还是慢衰落。

不同多径参数与信号经历的衰落类型之间的关系总结如图。

前面的一些理论的由来，以及公式推演，可以具体看看Smith Gorden的《Wireless Communication》

多径衰落信道统计特性

接下来，进入正题，我们的常见的衰落信道。接下来主要梳理的是瑞利衰落信道，因为莱斯信道只需要在瑞利信道的基础上增加一个直流分量就行（仅从公式来分析）。

瑞利衰落信道的一阶统计特性

当信道中不存在一个较强的直达径时，其信号包络服从是瑞利分布。

在移动无线信道中，Rayleigh分布是常见的用于描述平坦衰落信号或独立多径分量接收包络统计时变特性的一种分布类型。众所周知，两个正交的噪声信号之和的包络服从Rayleigh分布。Rayleigh分布的概率密度函数(pdf)为：

其中，$\sigma$是包络检波之前所接收的电压信号的rms值，${\sigma }^2$是包络检波之前的接收信号包络的时间平均功率。不超过某特定值R的接收信号的包络由相应的累积积分（CDF）给出：

Rayleigh分布的平均值r_mean为：

Rayleigh分布的方差为${\sigma }_r^2$，它表示信号包络的交流功率。表示为：

相位$\phi (t)$服从均匀分布的，即：

瑞利衰落信道的二阶统计特性

学习“概率论”时，我们便了解到往往一阶统计特性无法较好（准确）的来平均评价一个随机变量，当然随机过程也一样。

电平通过率（LCR）

电平通过率是指信道的幅度衰落电平在单位时间内以正斜率通过某规定电平R的次数的均值。每秒电平通过的数目为：

其中$\rho$

其中， r ̇是r(t)对时间的导数（即斜率）；$p(R,r)$是r=R处r与 r ̇ 的联合密度函数；$f_m$是最大多普勒频移；ρ是特定电平R相对于衰落包络的本地rms幅度进行归一化后的值。信号包络常常发生浅度衰落，偶尔也会发生深度衰落。

平均衰落持续时间（AFD）

平均衰落持续时间定义为接收信号低于某指定电平R的平均时间段的值，可以表示为：

由上式可知，AFD 与移动速率和载波频率成反比。可以通过计算平均衰落持续时间来近似确定在一次衰落期间最可能丢失的比特数。

以上的一些公式都是一些随机过程的一些数字特征，暂时有些不理解，其实没什么事，在后续的一些实践过程中，会慢慢的领会其中的含义。

所以看到现在有些灰心失意是正常哒，坚持看下去就是胜利呢；如果到现在感觉都没什么难度的话，那证明你的基本功很扎实呢，建议可以加快进度看下面的一些模型哈~

多径衰落信道模型建立

这里，我有查不少资料，有不少博主写的挺到位的，可参考
信道衰落随笔一，
瑞利信道：从原理到实现

小尺度平坦衰落信道中，由于移动体处于不断的运动状态，导致接收端接收到的是来自不同路径的多径信号。移动台的速率将会导致每条多径分量具有不同的多普勒频移，同时假设移动台接收的信号的入射角度在[ 0, 2π]间均匀分布。这样就产生了具有特定环境下的多普勒功率谱密度特性。Clarke信道模型描述的小尺度平坦衰落信道中，移动台接收信号具有基于散射方式的场强统计特性，这正符合市区移动通信环境不存在直射路径的特点,其包络统计特性服从瑞利分布，因此又称为瑞利衰落信道模型。一般情况下，在 瑞利衰落的状态下，多普勒功率谱具有Jakes 功率谱密度函数,即U型谱。

一般我们常见的衰落信道大多是基于多个不相干的高斯随机过程/有色高斯随机过程来构建。其中，瑞利信道和莱斯信道是由两组有色高斯来组成。

一般产生有色高斯噪声的方法有两类：第一类方法是成形滤波器法；第二类是正弦波叠加法。

基于clarke模型的成型滤波器法

成形滤波器法是将高斯白噪声输入成形滤波器来产生指定形状的多普勒功率谱，实部和虚部均是具有特定形状多普勒功率谱的高斯色噪声，因此包络服从瑞利分布。

有许多种具有可变增益和时延的Rayleigh 衰落仿真器，可以用它们来产生频率选择性衰落效应，如下图。

其中，下图整个频率选择性衰落效应的框图，可以理解成一个巨大的系统，该系统输入信号会收到三个变量的影响，分别是多径信号的相对时延扩展 以及各路径的增益 和多普勒扩展 （描述移动台的相对运动情况）。

基于Jakes模型的正弦波叠加法

正弦波叠加法利用确定性过程来模拟随机过程。根据中心极限定理，多个不同频率的正弦波叠加，幅度统计是服从高斯分布的。如果实部和虚部均用多个不同频率的正弦波叠加，则实部虚部均分别服从高斯分布。又根据瑞利分布的定义，两个高斯变量的平方和的根服从瑞利分布，因此包络服从瑞利分布。该方法易于实现，占用资源少，但是产生的多普勒功率谱是由多个离散频率点上的冲激构成的，并且相关性能不够理想。正弦波叠加法有很多种产生各正弦波相位、幅度参数的方法，比如等面积法，等距离法，蒙特卡罗法等。

式中

相移${\theta }_{(}i,n)$是$[0,2\pi ]$内均匀分布的随机变量；当$N_i\to \infty$时，$\Delta f_i\to 0$，这样就使频率成为连续分布的。

具体框图如下图

Jakes提出了基于正弦波叠加法的Jakes仿真模型。Jakes仿真模型是一种确定模型，产生的信号非广义平稳且不具各态历经性，其二阶统计特性与Clarke参考模型也相差较大。Jakes仿真模型虽然对Clarke参考模型实现了简化，降低了复杂度，提高了仿真效率，但却引入了广义非平稳性。其主要原因就在于Jakes对模型中的随机相移进行了确定化，同时造成了随机相移之间存在相关性。

下图是基于Jakes模型的正弦波叠加法，在工程中利用较多。

以类似的方式，可以将以上的模型封装为一个模块，用于多径衰落信道的模型建立。

多径衰落信道模型的仿真

接下来进行仿真了，其实这个仿真真心画了很久。一开始没怎么搞懂原理，就开始乱看学长学姐的代码，死磕了好久；后来意识到理论问题的时候，连忙去梳理这些知识，虽然距离仿真还是有一段距离，但是已经开始能搞懂大致的一些步骤，并且自己编写了基于clarke的成型滤波器算法。这阶段找了不少csdn的资料，也问了不少同学，都挺参差的。所以，不知不觉间，我就把以上的两个模型都建立和仿真了一下。

基于clarke的成型滤波器法仿真

首先先来看看基于clarke的成型滤波器法。

设定采样频率和多普勒频移扩展（这里主要是为了模拟一个时变的系统）。

这里构造了一个多普勒滤波器，不难发现复高斯噪声已经呈现“碗”型。下图更加明显一点（将其频域采样点数与输入信号统一化了）

然后，我们将以上的两组信号，分别ifft变回时域信号，进行平方求和后开根，得到瑞利衰落信道时域信号。时域波形将其与rms归一化后计算幅度增益，利用dB性质表示

绘制一节统计特性PDF，可能是因为点数不多，稍微有些畸变，整体服从瑞利分布趋势

二阶统计特性如下

可以发现这的电平通过率整体折半了，导致平均衰落时间是理论的两倍，我查了几次代码，感觉都查不出来，虽然LCR乘以2后波形呈现的很好，但是我总感觉是错的（理论支撑我写代码的，还是不能为了理论数据而造假，那就先放着等什么时候找到问题了，下面备注吧）

虽然这里的信道与往常的信道有些不同，这里的信道有些像我们通俗的乘性噪声（相对的，加性高斯白噪声，我们很熟悉，这是加性的，直接与信号相加就行），这里的信道需要与我们的输入信号（发射的信号），进行相乘；也可以通俗的理解为对输入信号的包络进行了变化。

根据卷积定理，时域相乘，频域卷积。所以我们需要知道产生的瑞利衰落信道的频谱图，这样方便我们来后续分析。以下是功率密度谱

看到这个波形的时候，我整个人惊讶到了。这个和我想象的频谱完全不同，怎么回事脉冲呢？不清楚你是不是和我有一样的疑惑。因为当时我是先做的Jakes模型的，在Jakes模型中其呈现的是碗的形状，而这里竟然是脉冲！

经过理论分析，我想这里是尤其合理性的。因为先前我们已经发现，这两条的任意一条之路的频谱都是碗的形状（因为被滤波了）。经过时域平方，相当于频域卷积，由于各频谱模的平方和，这必然是相当大的，自然呈现出脉冲状。

接下来，我们开始构建多径信道。假设多径个数为Num=50个，到来的每个信号的功率呈指数衰落，并且总功率和为1。我将多普勒频移最大带宽设置的较小，比如50Hz（假设其几乎静止），这样可以保证慢衰落（忘记概念的uu可以往前翻翻）。设置输入信号为一些cos信号的叠加（信号频率需要高于相干带宽118Hz）。

设置延时间隔（使相关带宽Bc=1422Hz大于输入信号之间的带宽），设置信道在平坦衰落处。以下是频域波形

平坦衰落的时域波形，基本没什么变化，就是增益变了

增大延时间隔（使相关带宽小于输入信号之间的带宽），设置信道在频率选择性衰落处。以下是频域波形

频率选择性衰落时域波形明显存在畸变

现在我们，保持平坦衰落，改变频域色散参数——多普勒频移最大带宽fm，原先是设置在50Hz，呈现出慢衰落，这时的时域波形如下

现在来看看快衰落的情况。

可以发现这里有畸变，但是畸变的非常微小。这个疑惑算是萦绕了我很久。我的理解是，由于之前瑞利信道的产生，其频谱的特性受多普勒频移扩展fm的影响太小了，所以对时域波形也影响较小。这里的思考有可能有误！（有大佬可以解答的话，希望得到您的指教）

基于Jakes的正弦波叠加法仿真

设定相同的参数，设置fm=50Hz。观察概率密度分布

由于fm设置太小，所以产生的正弦波叠加的个数较少，所以呈现出的分布有些畸形
当将fm调大后，PDF的波形明显更加逼近瑞利分布

此时的瑞利信道的PSD如下图

其二阶统计特性如下图，发现基本吻合理论值的情况。平均衰落时间在$\rho$较大的时候有偏差，这主要原因可能是测试时间较短所引起。

接下来，设置类似参数，首先设置在慢衰落，设置fm=5Hz(近乎理解为静止)。分析平坦衰落和频选衰落。
以下是平坦衰落

时域波形如下

当在频选衰落时

接下来，我设置在平坦衰落下，调节fm，呈现快衰落

此时时域波形明显发生畸变。从这里看出该模型的快慢衰落方面优于Clarke模型。

体会与心得

以上是全部的本次实验的仿真。

其实，还是有蛮多疑惑的地方的：

• 尤其是当快慢衰落与频选衰落纠缠在一起的时候。对于平坦，何时可以视为平坦，以及何时可以认为是快衰落（虽然有相关带宽和相干时间，但是对于代码编写者来说，无从得知其波形的正确性，只能凭借推理与想象了）
• 除此之外，就是频率精度对问题的复杂度的影响。因为Clarke模型的频谱近乎脉冲，所以影响较小；但是Jakes模型则会有影响，而且非常影响判断——比如明明处于平坦，但是却呈现出类似频选的特性，如下图

)

上图分别是频域和时域的图，根据相关带宽判断，应该是慢衰落、平坦衰落，但是其频谱却非常的奇怪；根本来说，是Jakes的频谱特性导致。

看到这里，如果你有什么想法或者建议的话，请在留言区留意，我们可以相互讨论。

代码

以下是主代码部分

clc;close all;clear;
%% 参数初始化
RayleighMode ='Jakes-Flat';%选择仿真模式'Forming-Filter'/'Jakes-Flat'
num = 50;%多径数目
power = 0:1/(num-1):1;
power = exp(-power);%各路径接收的功率
power = power/sum(power);
dt = 1e-6;%接收方采样间隔
delay = 1e-6*(0:(num-1));%各路径的附加时延
fsamp = 1/dt;%接收方采样频率
Last_time = 1e-1;%接收方信号采样持续时间
t = 0:dt:Last_time-dt;
n = length(t);
D = fsamp/n;%功率谱谱线间隔
fdmax = 50;%单径的多普勒扩展平移
tau = sum(power.*delay)/sum(power);%平均附加时延
tau2 = sum(power.*delay.^2)/sum(power);
sigma_tau = sqrt(tau2-tau^2);%rms时延扩展
Bc_09 = 1/50/sigma_tau;%频率相关函数为0.9时，相关带宽
Bc_05 = 1/5/sigma_tau;%频率相关函数为0.5时，相关带宽
Tc_05 = 0.423/fdmax;

%% 信号发生器
fs_limH = 3000;%有用信号基带高截止频率
fs_limL = 2500;%有用信号基带低截止频率
df = 10;%有用信号频率谱线间隔
fs = fs_limL:df:fs_limH;%有用信号频率
txSig = 0;
for ii = 1:length(fs)
    txSig = cos(2*pi*fs(ii)*t)+txSig;
end
%txSig = square(2*pi*fs*t);%square信号
% txSig = [1,zeros(1,n-1)];%脉冲信号
% flim_low = 0;
% flim_high = 100;
% txSig = FreSquare(n,fsamp,flim_low,flim_high);%频域矩形
% txSig = txSig.*sin(2*pi*fc*t);%sin信号

% txSig = ones(1,n);%全1信道

txSig = txSig';

%% 基于平坦瑞利信道的多径信号
output_signal=Frequency_selective_fading_Channel...
    (power,floor(delay/dt),txSig,fsamp,fdmax,RayleighMode);

%% 输入输出信号绘图
figure;hold on;box on;grid on
subplot(3,1,1);plot(t,txSig,'linewidth',2);title("输入信号");
subplot(3,2,5);plot(t,txSig,'linewidth',2);title("输入信号(放大)");xlim([0.0045,0.0055]);
subplot(3,1,2);plot(t,real(output_signal/(max(output_signal)-min(output_signal)))*2,'linewidth',2);title("频率选择性衰落后信号")
subplot(3,2,6);plot(t,real(output_signal/(max(output_signal)-min(output_signal)))*2,'linewidth',2);title("频率选择性衰落后信号(放大)");xlim([0.0045,0.0055]);
xlabel("时间/s");ylabel("幅度");

%% 输入输出信号的PSD
figure;box on;grid on;hold on;box on
fdata = fftshift(fft(txSig));
k=floor(-(n-1)/2: (n-1)/2);
D = fsamp/n;
subplot(3,1,1);plot(k*D,abs(fdata).^2,'linewidth',2);
xlim([-fs_limH*2,fs_limH*2]);
title("输入信号");
xlabel("频率/Hz");ylabel("PSD");
subplot(3,2,5);plot(k*D,abs(fdata).^2,'linewidth',0.5);
xlim([fs_limL*3/2-fs_limH/2,fs_limH*3/2-fs_limL/2]);
title("输入信号(正谱线)");
xlabel("频率/Hz");ylabel("PSD");

fout = fftshift(fft(output_signal));
k=floor(-(n-1)/2: (n-1)/2);
D = fsamp/n;
subplot(3,1,2);plot(k*D,abs(fout).^2,'linewidth',2);
xlim([-fs_limH*2,fs_limH*2]);
title("频率选择性衰落后信号");
xlabel("频率/Hz");ylabel("PSD");
subplot(3,2,6);plot(k*D,abs(fout).^2,'linewidth',0.5);
xlim([fs_limL*3/2-fs_limH/2,fs_limH*3/2-fs_limL/2]);
title("频率选择性衰落后信号(正谱线)");
xlabel("频率/Hz");ylabel("PSD");

%% 绘制外包络
output_signal_dB = 10*log10(output_signal/rms(output_signal));
figure;box on;grid on;
subplot(2,1,1);plot(t,output_signal_dB);
title("外包络衰落情况")
xlabel("Time/s");ylabel("Amplitude/dB");
subplot(2,1,2);plot(t,output_signal_dB);xlim([Last_time/2,Last_time/2+0.003]);
title("外包络衰落情况(x轴放大)")
xlabel("Time/s");ylabel("Amplitude/dB");

%% 二阶统计特性(计算瑞利分布的电平通过率（LCR）和平均衰落持续时间（AFD））
R_rms = sqrt(sum(abs(output_signal).^2)/n);
R_real = abs(output_signal)/R_rms;
R = 0:0.10:2;
rho = R;
rho2 = rho.^2;
N_R_theory = sqrt(2*pi)*fdmax*rho.*exp(-rho2);
N_R_real = zeros(1,length(rho));
ave_tau_theory = (exp(rho2)-1)./rho/fdmax/sqrt(2*pi);
ave_tau_real = zeros(1,length(rho));
for ii = 1:length(rho)
   for jj = 1:length(R_real)
       if(jj>=2)
           if((R_real(jj-1)<=R(ii))&&(R_real(jj)>=R(ii)))
               N_R_real(ii) = N_R_real(ii)+1;
           end
       end
       if(R_real(jj)<R(ii))
           ave_tau_real(ii) = ave_tau_real(ii)+1;
       end
   end
end
N_R_real = N_R_real/Last_time;
ave_tau_real = ave_tau_real/n./N_R_real;
figure;
subplot(2,1,1);box on;hold on;grid on;
plot(rho,N_R_theory,'bo-');
plot(rho,N_R_real,'r^-');
xlabel("\rho");ylabel("LCR");
legend("理论值","实际值");
subplot(2,1,2);box on;hold on;grid on;
plot(rho,ave_tau_theory,'bo-');
plot(rho,ave_tau_real,'r^-');
xlabel("\rho");ylabel("AFD");
legend("理论值","实际值");

下面的是主代码中用到的一些函数

function output_signal=Frequency_selective_fading_Channel(power,delay,txSig,fsamp,fdmax,RayleighMode)
%power PDP对应的power 单位W 
%delay PDP对应的delay 单位us 因此此处采样频率需为1MHZ 
%txSig 信道输入信号 输入为列向量 
%fsamp 输入信号的抽样频率 
fs=fsamp; %本程序的抽样频率 
txSig_len=length(txSig); 
ratio=ceil(fs/fsamp); %过采样率
col=ratio*txSig_len; %仿真矩阵列数 
row=length(delay); %仿真矩阵行数
txSig_sim=zeros(row,col);%输入信号模拟
h_sim=zeros(row,col); %信道模拟
sigma0_2=1/2; %噪声功率
T_interval=1/fs; %采样时间间隔

txSig_spread=ones(1,ratio).*txSig(1); 
for i=2:txSig_len 
    txSig_spread=[txSig_spread,ones(1,ratio).*txSig(i)]; 
end
txSig_sim=repmat(txSig_spread,row,1); 
 
switch(RayleighMode(1:2))
    case 'Ja'
        N0=500;%正弦波叠加数，一般大于50
        for ii = 1:row
            h_sim(ii,:)=Rayleigh_Channel(fdmax,sigma0_2,N0,T_interval,1,col); %基于'Jakes-Flat'
        end
    case 'Fo'
        for ii = 1:row
            h_sim(ii,:)=Rayleigh_Channel1(fdmax,sigma0_2,T_interval,1,col); %基于'Forming-Filter'
        end
    otherwise
        error('Unknown Rayleigh  Simulation Model in Frequency_selective_fading_Channel()');
end

%% 对信号按照延迟进行移位 
figure;hold on% 输入信号和信道的PSD
n = length(txSig);
k=floor(-(n-1)/2: (n-1)/2);
fdata = fftshift(fft(txSig));
D = fsamp/n;
subplot(2,1,1);plot(k*D,abs(fdata).^2,'linewidth',2);

D = fs/n;
subplot(2,1,2);hold on;
for i=1:row 
    txSig_sim(i,:)=[zeros(1,delay(i)),txSig_sim(i,1:(end-delay(i)))]; %信号时延
    resul_mat(i,:) = txSig_sim(i,:).*h_sim(i,:)*power(i);
    fh = fftshift(fft(h_sim(i,:)));
    plot(k*D,abs(fh).^2,'linewidth',2);
end
xlabel("Frequency/Hz");ylabel("PSD");
% legend("输入信号","瑞利信道")

%%输出信号
output_signal_spread=sum(resul_mat,1); %多径叠加
index=(1:txSig_len)*ratio; %符号同步
output_signal=output_signal_spread(index);%下采样 
end

function T_out=Rayleigh_Channel(fdmax,sigma_0_2,N0,T_interval,row,col) 
%fdmax=40; 指定多普勒扩展 
%N0=25; 指定正弦波叠加的数目 
%% 此模型基于《移动传播环境》P215 jakes模型 以及P220模型的改进 

%% 初始化各参数
t=T_interval.*[1:(row*col)]; %仿真时长与间隔 
N=4*N0+2; %总的仿真点数
n=1:N0; 
wd =2*pi*fdmax;
wn = wd.*cos(2*pi/N*n);
beta_n=pi.*(n/(N0+1)); 
phase_n=2*pi*rand(1,N0); 
Tc=zeros(size(t)); 
Ts=zeros(size(t)); 
for i=1:N0 
    Tc = Tc+cos(beta_n(i))*cos(wn(i)*t+phase_n(i));
    Ts = Ts+sin(beta_n(i))*cos(wn(i)*t+phase_n(i));
end
Tc = Tc*2+sqrt(2)*cos(wd*t);
Ts = Ts*2;
T_out=sqrt(sigma_0_2)*reshape((Tc+1j*Ts),row,col)/sqrt(2*N0+1); 

%% CDF/PDF
% figure;box on;
% ksdensity(abs(Tc+1*j*Ts));%PDF
% xlabel("x");ylabel("PDF");

%% PSD
% n = length(T_out);
% k=floor(-(n-1)/2: (n-1)/2);
% fdata = fftshift(fft(T_out));
% fsamp = 1e6;
% D = fsamp/n;
% figure;
% subplot(2,1,1);plot(k*D,abs(fdata).^2,'linewidth',2);
% ylabel("PSD");xlabel("频率");title("基于正弦波叠加的瑞利信道PSD")
% subplot(2,1,2);plot(k*D,abs(fdata).^2,'linewidth',2);xlim([-2*fdmax,2*fdmax])
% ylabel("PSD");xlabel("频率");title("基于正弦波叠加的瑞利信道PSD(x轴放大)")
end

function T_out=Rayleigh_Channel1(fdmax,sigma_0_2,T_interval,row,col) 
% clc;close all;clear;
%% 测试参数初始化
% fm=10; %指定多普勒扩展Hz
% N = 1000;%采样点数
% sigma_0_2 = 1/2;%噪声功率
% df = 2*fm/(N-1);%相邻谱线频率间隔Hz
% dT = 1/df;%衰落波形的时间周期s
% fs 采样频率

%% 函数参数初始化
fm = fdmax;%多普勒频移宽度-fm~fm
N = row*col;%总采样点数
dT = T_interval;%采样时间间隔
fs = 1/dT;%采样频率
N0 = floor(fdmax/(fs/N));%fd范围内的频点数
f = fdmax*2*floor(-(N0-1)/2: (N0-1)/2)/N0;
index = find(f==fm);%找到多普勒边界值
f(index)=f(index)-(fs/N)/100;%避免无穷大的异常值
index = find(f==-fm);%找到多普勒边界值
f(index)=f(index)+(fs/N)/100;%避免无穷大的异常值

%% 复高斯噪声发生器(从频域产生)
Sig1 = randn(1,floor(N0/2))+i*randn(1,floor(N0/2));%正频域的复高斯变量
Sig2 = randn(1,floor(N0/2))+i*randn(1,floor(N0/2));%正频域的复高斯变量
if(mod(N0,2)==1)%需要补零频
    Sig1 = [conj(Sig1(end:-1:1)),0,Sig1];%正负频域的复高斯变量(直流分量设为0)
    Sig2 = [conj(Sig2(end:-1:1)),0,Sig2];%正负频域的复高斯变量(直流分量设为0)
else
    Sig1 = [conj(Sig1(end:-1:1)),Sig1];%正负频域的复高斯变量(直流分量设为0)
    Sig2 = [conj(Sig2(end:-1:1)),Sig2];%正负频域的复高斯变量(直流分量设为0)
end

%% 基带多普勒滤波器(-fm/2~fm/2之间)
SF1 = Sig1.*sqrt(1.5/pi/fm./sqrt(1-(f/fm).^2));
SF2 = Sig2.*sqrt(1.5/pi/fm./sqrt(1-(f/fm).^2));

% figure;box on;grid on;hold on;
% set (gca,'unit','normalized','position',[0.1,0.3,0.8,0.4] );
% subplot(2,1,1);%stem(f,real(Sig1),'linewidth',2);title("复高斯线谱")
% hold on;plot(f,sqrt(1.5/pi/(fm/2)./sqrt(1-(f/(fm/2)).^2)));
% subplot(2,1,2);stem(f,abs(SF1),'linewidth',2);title("多普勒滤波后线谱")

%% IFFT
SF1_Shift = fftshift(SF1);
SF1_Expand = [SF1_Shift(1:floor(N0/2)),zeros(1,N-N0),SF1_Shift(floor(N0/2)+1:end)];
SF2_Shift = fftshift(SF2);
SF2_Expand = [SF2_Shift(1:floor(N0/2)),zeros(1,N-N0),SF2_Shift(floor(N0/2)+1:end)];
sf1 = ifft(SF1_Expand);
sf2 = ifft(SF2_Expand)*exp(-1i*pi/2);%正交分量相移90度

% figure;
% k=floor(-(N-1)/2: (N-1)/2);
% D = fs/N;
% subplot(2,1,1);plot(k*D,abs(fftshift(SF1_Expand)),'linewidth',2);xlabel("频率/Hz");
% title("经过多普勒滤波器后的谱线")
% subplot(2,1,2);plot(k*D,abs(fftshift(SF1_Expand)),'linewidth',2);xlabel("频率/Hz");
% xlim([-2*fm,2*fm]);
% title("经过多普勒滤波器后的谱线(放大)")

%% 瑞利仿真衰落信号
sf = sqrt(real(sf1).^2+real(sf2).^2);
sf = sf*sqrt(sigma_0_2/mean(sf.*conj(sf)));%功率调节
% figure;box on
% sf_dB = 10*log10(sf/rms(sf));
% plot([1:N]*dT,sf_dB);xlabel("时间/s");ylabel("归一化增益/dB")
% title("瑞利衰落包络信号")

%% CDF/PDF
% figure;box on;
% ksdensity(real(sf));%PDF
% xlabel("x");ylabel("PDF");

%% PSD
% figure;box on;grid on;hold on;box on
% fdata = fftshift(fft(sf));
% k=floor(-(N-1)/2: (N-1)/2);
% D = fs/N;
% subplot(2,1,1);plot(k*D,abs(fdata).^2,'linewidth',2);ylabel("PSD");xlabel("频率/Hz");title("瑞利衰落信号PSD")
% subplot(2,1,2);stem(k*D,abs(fdata).^2,'linewidth',2);ylabel("PSD");xlabel("频率/Hz");title("瑞利衰落信号PSD(x轴放大)");
% xlim([-fs/N*5,fs/N*5]);

%% row条信道
T_out=reshape(sf,row,col);