滤波白噪声模型

产生滤波白噪声信道模型，画出信道幅度系数及其包络的概率密度函数与相位的概率密度函数

clear,clf
fm=100;%最大多普勒频率
scale=1e-6;%微秒量级
ts_mu=50;ts=ts_mu*scale;fs=1/ts;%采样时间与采样频率
Nd=1e6;%采样个数
%获得复信道系数
[h,Nfft,Nifft,doppler_coeff]=FWGN_model(fm,fs,Nd);
subplot(211)
plot([1:Nd]*ts,10*log10(abs(h)))%画信道幅度系数
str=sprintf('Clarke/Gan Model,f_m=%d[Hz],T_s=%d[us]',fm,ts_mu);
title(str),axis([0 0.5 -30 5])
subplot(223)
histogram(abs(h),50)%画信道包络概率密度函数，包络一维概率密度函数呈瑞利分布
subplot(224)
histogram(angle(h),50)%画信道相位概率密度函数，相位一维概率密度函数呈均匀分布

%滤波白噪声信道子程序（Clarke/Gan模型），子程序名称：FWGN_model.m
function [h,Nfft,Nifft,doppler_coeff] = FWGN_model(fm,fs,N)
%输入参数：fm为最大多普勒频率，fs为采样频率，N为采样个数
%输出参数：h为复信道系数
Nfft =2^nextpow2(2*fm/fs*N);% nextpow2是用来求指数的，这个指数（假定为a）需满足两个条件：（1）2的a次方大于等于输入值；（2）2的（a-1）次方小于输入值。
Nifft = ceil(Nfft*fs/(2*fm));
%产生独立的复高斯随机过程，参见图2-12
GI=randn(1,Nfft);GQ=randn(1,Nfft);
%求实信号的FFT，以获得Hermitian对称
CGI=fft(GI);CGQ=fft(GQ);%将高斯随机过程转换到频域
%多普勒谱产生，仿真多普勒滤波器
doppler_coeff=Doppler_spectrum(fm,Nfft);
%将转换到频域的高斯过程加入到多普勒滤波器，频域内乘积
f_CGI=CGI.*sqrt(doppler_coeff);
f_CGQ=CGQ.*sqrt(doppler_coeff);
%补零，使多普勒滤波器输出数据长度为Nifft，以备求ifft
Filtered_CGI=[f_CGI(1:Nfft/2) zeros(1,Nifft-Nfft) f_CGI(Nfft/2+1:Nfft)];
Filtered_CGQ=[f_CGQ(1:Nfft/2) zeros(1,Nifft-Nfft) f_CGQ(Nfft/2+1:Nfft)];
%求ifft，将多普勒滤波器输出频域信号转换到时域
hI=ifft(Filtered_CGI);hQ=(Filtered_CGQ);
%计算实部的平方加虚部的平方开平方，即呈瑞利分布的包络
rayEnvelope=sqrt(abs(hI).^2+abs(hQ).^2);
%计算包络的方均根值
rayRMS=sqrt(mean(rayEnvelope(1:N).*rayEnvelope(1:N)));
%图2-12最终输出的Clarke/Gan模型信道冲激响应
h=complex(real(hI(1:N)),-real(hQ(1:N)))/rayRMS;

%多普勒子程序，子程序名称：Doppler_spectrum.m
%对经典多普勒谱的仿真，参见式（2-25）
function y=Doppler_spectrum(fd,Nfft)
%输入参数：fd为最大多普勒频移，Nfft为频域样值点个数
%输出参数：y返回多普勒谱
df=2*fd/Nfft;%计算频率间隔
%计算f=0时多普勒谱
f(1)=0;y(1)=1.5/(pi*fd);
%计算其他频率多普勒谱。计算多普勒谱从第2个样值到Nfft/2个样值以及第Nfft/2+2个样值到第Nfft个样值
for i = 2:Nfft/2
    f(i)=(i-1)*df;%根据频率序号及频率间隔计算频率
    y([i Nfft-i+2])=1.5/(pi*fd*sqrt(1-(f(i)/fd)^2));
    %计算式(2-25),注意randn产生的高斯噪声方差为1
end
%计算多普勒谱的Nfft/2+1点样值。用四点样值构成多项式，之后求出多项式在Nfft/2+1点的样值。
nFitPoints=3;
kk=[Nfft/2-nFitPoints:Nfft/2];%四点样值序号为Nfft/2-3到Nfft/2
polyFreq=polyfit(f(kk),y(kk),3);%构成多项式
y((Nfft/2)+1)=polyval(polyFreq,f(Nfft/2)+df);%从多项式求Nfft/2点值

在我电脑上的matlab中遇到了这个警告，我将阶数修改为2虽然没有了警告，但图形和参考书上的还是对不上，目前还没找到那里的问题…

参考《MIMO-OFDM系统原理、应用及仿真》一书