%main.m
clc;
LengthOfSignal=10000; %信号长度
fm=500;    %最大多普勒频移？相关文献应该有估算公式
fc=5000;   %信道载波频率
t=1:LengthOfSignal;  % SignalInput=sin(t/100);
%DSB调制
SignalInput=sin(t/50);%+cos(t/65);  %调制信号
c=cos(0.2*pi*t);%载波信号
y_in=SignalInput.*c;%调制
delay=[0 31 71 109 173 251];%10ns
power=[0 -1 -9 -10 -15 -20];   %dB
y_in=[zeros(1,delay(6)) y_in];   %为时移补零
y_out=zeros(1,LengthOfSignal);   %存放经信道未解调的信号（现为无输入信号
%时的输出信号）
%y_out_end最终解调后信号
%多路径衰落
for i=1:6%图4
    f=1:2*fm-1;
    Rayl;
    y_out=y_out+r.*y_in(delay(6)+1-delay(i):(delay(6)+LengthOfSignal-delay(i)))*10^(power(i)/20);
end;
% S(t)*cos(w*t)=m(t)*cos(w*t)*cos(w*t)=0.5*m(t)*(1+cos(2*w*t))
%用一个低通滤波器将上式中的第一项和第二项分离，无失真的恢复出原始的调制信号。
%这种调制方法又称为同步解调或相干解调
%同步解调
y_out_end=y_out.*c;%同步解调或相干解调
%低通滤波
wp=0.1*pi;ws=0.12*pi;Rp=1;As=15;
[N,wn]=buttord(wp/pi,ws/pi,Rp,As);
[b,a]=butter(N,wn);
y_out_end =filter(b,a,y_out_end);%滤波
y_out_end =2* y_out_end;%恢复幅度
%原信号的频谱
K=fft(SignalInput);
%DSB调制后信号的频谱
L=fft(y_in);
%y_out的频谱（含包络）
M=fft(y_out);
%最终解调的频谱
N=fft(y_out_end);
%输出
figure(1);
subplot(4,2,1);
plot(SignalInput(delay(6)+1:LengthOfSignal));axis([0,3000,-2,2]);
title('原始输入信号');
subplot(4,2,2);
plot(abs(fftshift(K))); axis([4900,5100,0,6000]);
title('原始输入信号的频谱 ');
subplot(4,2,3);
plot(y_in(delay(6)+1:LengthOfSignal));axis([0,3000,-2,2]); %去除时延造成的空白信号
title(' 进入瑞利信道前，DSB调制后的信号');
subplot(4,2,4);
plot(abs(fftshift(L))); axis([3500,6500,0,3000]);
title('进入瑞利信道前，DSB调制后的信号的频谱 ');
subplot(4,2,5);
plot(y_out(delay(6)+1:LengthOfSignal));axis([0,3000,-0.08,0.08]); %去除时延造成的空白信号
title('经瑞利信道后，DSB解调前的信号');
subplot(4,2,6);
plot(abs(fftshift(M)));axis([3500,6500,0,100]);
title('经瑞利信道后，DSB解调前的信号的频谱');
subplot(4,2,7);
plot(y_out_end(delay(6)+1:LengthOfSignal));axis([0,3000,-0.08,0.08]);%去除时延造成的空白信号
title('最终解调后的信号');
subplot(4,2,8);
plot(abs(fftshift(N)));axis([4900,5100,0,200]);
title('最终解调后的信号的频谱 ');
figure(2);
subplot(3,1,1);
hist(r,256);%绘制直方图
title('瑞利信道的幅度分布')
subplot(3,1,2);
hist(angle(r0));
title('瑞利信道的相位分布');
subplot(3,1,3);
plot(Sf1);
title('多普勒滤波器的频响特性');

3.2瑞利衰落信道的matlab实现
%Rayl.m 参考【1】
f=1:2*fm-1; %通频带长度
y=1.5./((1-((f-fm)/fm).^2).^(1/2))/pi/fm; %多普勒功率谱(基带)图3
Sf=zeros(1,LengthOfSignal);
Sf1=y;
Sf(fc-fm+1:fc+fm-1)=y;%(把基带映射到载波频率)
x1=randn(1,LengthOfSignal);
x2=randn(1,LengthOfSignal);
nc=ifft(fft(x1+1i*x2).*sqrt(Sf)); %同相分量
%首先产生独立的复高斯噪声的样本，并经过FFT后形成频域的样本，
%然后与S（f）开方后的值相乘，以获得满足多普勒频谱特性要求的信号，
%经IFFT后变换成时域波形，再经过平方，将两路的信号相加并进行开方运算后，形成瑞利衰落的信号r(t)
x3=randn(1,LengthOfSignal);
x4=randn(1,LengthOfSignal);
ns=ifft(fft(x3+1i*x4).*sqrt(Sf)); %正交分量
r0=(real(nc)+1i*real(ns)); %瑞利信号
r=abs(r0); %瑞利信号幅值（nc、ns分别为窄带高斯过程的同相和正交支路的基带信号）