matlab实现频域瑞利(Rayleigh)信道仿真

clear;

fc=input('请输入载波频率(Hz)：');

v=input('请输入接收机速度(Km/h)：');

v1=v*1000/3600;%接收机速率[km/h]

c=300*10^6;%光速

fm=fc*(v1/c);

N =input('请输入频域点数：');

% 产生多普勒功率谱

gap = 2*fm/(N-1);

T = 1/gap;

sf0 = 1.5/(pi*fm);

for n = 1:(N-2)/2

sf(n) = 1.5/(pi*fm*sqrt(1-(n*gap/fm)^2));

end

SEf = [fliplr(sf),sf0,sf]; %合成全频段的多普勒功率谱

figure(1);

plot(SEf);

title('多普勒功率谱');

xlabel('f');

ylabel('size');

grid;

% 产生两个正态分布噪声源

Gauss_time1=randn(1,N-1);

Gauss_time2=randn(1,N-1);

GaussN1=fft(Gauss_time1);

GaussN2=fft(Gauss_time2);

% 产生瑞利衰落信道

x = ifft(sqrt(SEf).*GaussN1);

y = ifft(sqrt(SEf).*GaussN2);

rayleigh_amp = sqrt(abs(x).^2+abs(y).^2);

rayleigh_db = 20*log10(rayleigh_amp); %用dB表示瑞利信号

figure(2);

plot(rayleigh_db);

% axis([0 140 -100 20]);

title('瑞利信号衰落');

xlabel('N');

ylabel('dB');

grid;

figure(3)

r = sqrt(0.5*(real(Gauss_time1).^2 + real(Gauss_time2).^2));

step = 0.1; range = 0:step:3;

h = hist(r, range);

fr_approx = h/(step*sum(h));

fr = (range/0.5).*exp(-range.^2);

plot(range, fr_approx,'r', range, fr,'k');

title('瑞利概率密度');

xlabel('rms');

ylabel('p(r)');

grid;

figure(4)

subplot(2,2,1)

plot(Gauss_time1);

title('时域高斯信号1');

xlabel('N');

ylabel('V');

grid;

subplot(2,2,2)

plot(Gauss_time2);

title('时域高斯信号2');

xlabel('N');

ylabel('V');

grid;

subplot(2,2,3)

plot(GaussN1);

title('频域复数高斯信号1');

xlabel('实部');

ylabel('虚部');

grid;

subplot(2,2,4)

plot(GaussN2);

title('频域复数高斯信号2');

xlabel('实部');

ylabel('虚部');

grid;

figure(5)

subplot(2,1,1)

plot(sqrt(SEf).*GaussN1);

title('高斯噪声与多普勒功率谱相乘1');

xlabel('实部');

ylabel('虚部');

grid;

subplot(2,1,2)

plot(sqrt(SEf).*GaussN2);

title('高斯噪声与多普勒功率谱相乘2');

xlabel('实部');

ylabel('虚部');

grid;

figure(6)

subplot(2,1,1)

plot(x);

title('IFFT1');

xlabel('实部');

ylabel('虚部');

grid;

subplot(2,1,2)

plot(y);

title('IFFT1');

xlabel('实部');

ylabel('虚部');

grid;

figure(7)

subplot(2,1,1)

plot(abs(x).^2);

title('求模平方1');

xlabel('N');

ylabel('V');

grid;

subplot(2,1,2)

plot(abs(y).^2);

title('求模平方2');

xlabel('N');

ylabel('V');

grid;

figure(8)

subplot(2,1,1)

plot(abs(x).^2+abs(y).^2);

title('平方后相加');

xlabel('N');

ylabel('V');

grid;

subplot(2,1,2)

plot(sqrt(abs(x).^2+abs(y).^2));

title('相加后开方');

xlabel('N');

ylabel('V');

grid;

%求均值、方差、均方差

sum1=0;

sum2=0;

sum3=0;

for n=1:length(rayleigh_amp)

sum1=sum1+rayleigh_amp(n);

end

ave=sum1/length(rayleigh_amp) %求瑞利衰减均值

for n=1:length(rayleigh_amp)

sum2=sum2+rayleigh_amp(n)^2;

sum3=sum3+(rayleigh_amp(n)-ave)^2;

end

z=sum2/length(rayleigh_amp);

fc=sum3/lengt

h(rayleigh_amp) %求出瑞利衰减方差

rms=sqrt(z) %求出瑞利衰减均方差

%求电平通过率、平均衰落持续时间和误比特率

i=0;

j=0;

e=0;

q=0;

for n=1:length(rayleigh_amp)-1

if(((rayleigh_amp(n+1)-rayleigh_amp(n))>=0)&(rayleigh_amp(n+1)>=rms)&(rayleigh_amp(n)<=rms))

i=i+1;

end

end

LCR=i/T %求出电平通过率

for m=1:length(fr_approx)

if(rayleigh_amp(m)<=rms)

j=j+1;

e=fr_approx(m)+e;

end

end

time=j*T/length(fr_approx);

AFD=e/LCR %求出平均衰落持续时间

bps=length(rayleigh_amp)/T %计算数据速率

BER=LCR/bps %求出误比特率