返回信息流从网上找到的一个基于clark模型的瑞利信道的程序：
function r = Channel_rayleigh( fd, fs, Ns )
% r = rayleigh(fd,fs,N)
%
% A Rayleigh fading simulator based on Clarke's Model
% Creates a Rayleigh random process with PSD determined
% by the vehicle's speed.
%
%INPUTS:
% fd = doppler frequency
% set fd = v*cos(theta)/lambda
% v = velocity (meters per second)
% lambda = carrier wavelength (meters)
% theta = angle w/ respect to tangent (radians).
% fs = sample frequency (Samples per second)
% Ns = number of samples of the Rayleigh fading
% process to produce
%
% OUTPUTS:
% r = row vector containing Ns samples of the Rayleigh
% fading process
%
%v = input('enter the velocity (meters per second): ')
%lambda = input('enter the carrier wavelength (meters): ')
%theta = input('enter the angle w/ respect to tangent (radians): ')
%fd = v*cos(theta)/lambda
%fd = input('doppler frequency: ');
%fs = input('enter the sample frequency (Samples per second): ');
%Ns = input(' enter the number of samples of the Rayleigh fading process to produce: ');
N = 8;
while(N)
if (N < 2*fd*Ns/fs)
N = 2*N;
else
break;
end
end % number of ifft points (for smoothing)
N_inv = ceil(N*fs/(2*fd));
% determine the frequency spacings of signal after FFT
delta_f = 2*fd/N;
% determine time spacing of output samples
delta_T_inv = 1/fs;
%%%%%%%%%%% Begin Random Input Generation %%%%%%%%%%%%
% fprintf( 'Generating Input\n');
% generate a pair of TIME DOMAIN gaussian i.i.d. r.v.'s
I_input_time1 = randn(1,N);
Q_input_time1 = randn(1,N);
% take FFT
I_input_freq = fft(I_input_time1);
Q_input_freq = fft(Q_input_time1);
%%% Generate Doppler Filter's Frequency Response %%%
% fprintf( 'Generating Doppler Filter Function\n');
% filter's DC component
SEZ(1) = 1.5/(pi*fd);
% 0 < f < fd
for j=2:N/2
f(j) = (j-1)*delta_f;
SEZ(j) = 1.5/(pi*fd*sqrt(1-(f(j)/fd)^2));
SEZ(N-j+2) = SEZ(j);
end
% use a polynomial fit to get the component at f = fd
% p = polyfit( f(N/2-3:N/2), SEZ(N/2-3:N/2), 3);
% SEZ(N/2+1) = polyval( p, f(N/2)+delta_f );
%k = N/2 - 1;
%For deep fading, picking tangential point in filter
k = 3;
p = polyfit( f(N/2-k:N/2), SEZ(N/2-k:N/2), k );
SEZ(N/2+1) = polyval( p, f(N/2)+delta_f );
%%%%%%%% Perform Filtering Operation %%%%%%%%%%%%%%
% fprintf( 'Computing Output\n' );
% pass the input freq. components through the filter
I_output_freq = I_input_freq .* sqrt(SEZ);
Q_output_freq = Q_input_freq .* sqrt(SEZ);
% take inverse FFT
I_temp = [I_output_freq(1:N/2) zeros(1,N_inv-N) I_output_freq(N/2+1:N)];
I_output_time =real(ifft(I_temp));
Q_temp = [Q_output_freq(1:N/2) zeros(1,N_inv-N) Q_output_freq(N/2+1:N)];
Q_output_time =real(ifft(Q_temp));
% take magnitude squared of each component and add together
for j=1:N_inv
r(j) = sqrt( (I_output_time(j))^2 + (Q_output_time(j))^2);
end
% normalize and compute rms level
rms = sqrt( mean( r.*r ) );
r=r(1:Ns)/rms;% r=(I_output_time+i*Q_output_time)/rms;
% r=r(1:Ns);
% s=real(r);
% t=imag(r);
% for ii=1:Ns
% r(ii)=sqrt((s(ii))^2 + (t(ii))^2);
% end
% r = r(1:Ns);
标出红色的部分不太明白，麻烦知道的指导一下，多谢了~~
这是一条镜像帖。来源：北邮人论坛 / communications / #4250同步于 2007/5/28
该镜像源已超过 30 天没有更新，可能在源站已被删除。
Communications机器人发帖
Channel_rayleigh函数的问题
kitty101103
2007/5/28镜像同步10 回复
订阅后，新回复会通过你的通知中心匿名送达。
9 条回复
N = 8;
while(N)
if (N < 2*fd*Ns/fs)
N = 2*N;
else
break;
end
end % number of ifft points (for smoothing)
这段相当于：N = 2^ceil(log2(2*fd*Ns/fs));
目的是得到不小于Doppler对应Nyquest Sampling Frequency的2的整数幂，作为Doppler频域采样点数
这种方法是采用频域滤波的方法生成衰落。
瑞利平坦衰落实际上是有色的高斯随机过程，根据Clarke的Density Scattering模型，这种衰落的包络的互相关函数是Modified Bessel function on zero order，对应的功率谱密度函数是Classic spectrum(U型谱)。
既然是有色的高斯过程，那么最直接的生成方法就是将白色的高斯过程通过一个滤波器，对这个高斯过程进行“着色”。
显然，这个用来着色的滤波器就是你希望的功率谱密度开根号，因为着色是对幅度操作。
你现在需要的是U型谱，因此：
SEZ(1) = 1.5/(pi*fd);
% 0 < f < fd
for j=2:N/2
f(j) = (j-1)*delta_f;
SEZ(j) = 1.5/(pi*fd*sqrt(1-(f(j)/fd)^2));
SEZ(N-j+2) = SEZ(j);
end
这一段是生成U型谱。
由于在Doppler处的功率谱密度趋于无穷大，因此：
k = 3;
p = polyfit( f(N/2-k:N/2), SEZ(N/2-k:N/2), k );
SEZ(N/2+1) = polyval( p, f(N/2)+delta_f );
这段对Doppler处进行多项式插值；
现在你就得到了滤波器的频域传递函数，此时只需将输入信号转到频域，然后相乘，在将结果转回时域，就是你希望的时域的衰落样点。
因此，此段
I_input_time1 = randn(1,N);
Q_input_time1 = randn(1,N);
% take FFT
I_input_freq = fft(I_input_time1);
Q_input_freq = fft(Q_input_time1);
将白高斯过程变到频域，然后：
I_output_freq = I_input_freq .* sqrt(SEZ);
Q_output_freq = Q_input_freq .* sqrt(SEZ);
进行相乘，然后：
I_temp = [I_output_freq(1:N/2) zeros(1,N_inv-N) I_output_freq(N/2+1:N)];
I_output_time =real(ifft(I_temp));
Q_temp = [Q_output_freq(1:N/2) zeros(1,N_inv-N) Q_output_freq(N/2+1:N)];
Q_output_time =real(ifft(Q_temp));
把结果再转回时域。
% take magnitude squared of each component and add together
for j=1:N_inv
r(j) = sqrt( (I_output_time(j))^2 + (Q_output_time(j))^2);
end
% normalize and compute rms level
rms = sqrt( mean( r.*r ) );
r=r(1:Ns)/rms;% r=(I_output_time+i*Q_output_time)/rms;
是为了将输出样点的功率进行归一化，因为衰落的功率是1。
附一篇信道估计文章：
附件(391.6KB) f_channel_estimation_techniques_in_multipath_fading_CDMA.pdf
谢谢你的热心解答，下面的这段：
SEZ(1) = 1.5/(pi*fd);
% 0 < f < fd
for j=2:N/2
f(j) = (j-1)*delta_f;
SEZ(j) = 1.5/(pi*fd*sqrt(1-(f(j)/fd)^2));
SEZ(N-j+2) = SEZ(j);
end
这一段是生成U型谱。
由于在Doppler处的功率谱密度趋于无穷大，因此：
k = 3;
p = polyfit( f(N/2-k:N/2), SEZ(N/2-k:N/2), k );
SEZ(N/2+1) = polyval( p, f(N/2)+delta_f );
这段对Doppler处进行多项式插值；
生成u型谱的程序和多项式插值部分还是不太明白，能再详细介绍下吗？麻烦了~
还有，cdma系统在瑞丽信道下是不是容易有平板效应（就是说不考虑信道估计的话，系统的snr-ber曲线当信噪比较大时，比较平稳）？
如果你不做信道估计的话，你就没办法做RAKE，也就是说，cdma最关键的技术你没办法利用，性能将很差。因为多径之间的干扰很大很大，而且，最关键的是，多径之间的干扰和信噪比没什么关系，因为它是信号之间的干扰，就算你提高信号的发射功率，它们之间的干扰还是很大。
此外，fd是Doppler，如果不引入它的话，怎么做Rayleigh信道呢？
哦，原来信道估计是非常重要的，我主要没怎么研究这方面，不了解，谢谢了~~
【 在 ismydoom 的大作中提到: 】
: 而fs是采样率。U型滤波器是用数字FIR实现的，你需要进行采样，fs就是采样率。
关于LS（Linear Smoothing）channel estimation
LS估计最简单的办法就是对导频信号进行低通滤波
设受到的导频信号为y(k) = x(k)+n(k),都是复信号；接收机信道估计：
y_est = sum(y(k)*x_local(k))/N；
其中，x_local是本地产生的导频信号，满足x_local(k)*x(k)=1
因此y_est=sum(1+n(k)*x_local(k))/N=1+sum(n(k)*x_local(k))/N；
设z=sum(n(k)*x_local(k))/N，由于x_local是PN序列，易知z~N(0, (sigma^2)/N)，
也就是说，通过低通滤波器（线性平滑）后，信噪比提高了N倍，然后用y_est作为信道的冲激响应，进一步进行RAKE的MRC等操作