matlab 实现频域瑞利(Ra yleigh)信道仿真 %接收机速率[km/h] %光速 clc; clear; fc=input('请输入载波频率(Hz)：'); v=input('请输入接收机速度(Km/h)：'); v1=v*1000/3600; 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/length(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 %求出误比特率