clarke_mod.m

代码为实现了信道CLARKE模型的仿真代码

function imp_rsp_c=clarke_mod(fs,Ns,Fdm)
%function imp_rsp_c=clarke_mod(Ns,Fdm)
%
%   this function produce the impulse respond of the channel based by
%       Clarke model.This Rayleigh fading simulator creates a Rayleigh 
%       random process with PSD determined by the vehicle's speed.
%   The algorithm is base on Smith simulating tool.
%
%   fs:  sample frequency (Samples per second)
%   Ns:  number of samples of the Rayleigh fading process to produce;
%   Fdm: maximum doppler frequency;
%       set Fm = v*cos(alpha)/lambda
%        v = velocity ,lambda = carrier wavelength ,alpha = angle w/ respect to tangent (radians).
%
%   imp_rsp_c:a row vector containing Ns samples of the Rayleigh fading process .

%Ns=100000;Fdm=1;

%fs=1000000;   %here we sets the sample frequency 1.
N=8;
while(N)
   if (N<2*Fdm*Ns/fs) 
      N = 2*N;
   else
      break;
   end
end 
%   1---the number of frequency domain points Ns used to represent sqrt(SEZ(f)) and the muximum Doppler freqnency
%       shift Fdm are input parameters.
%   2---computer the frequency spacing between adjacent spectral lines---delta_f

% number of ifft points (for smoothing)
N_inv = ceil(N*fs/(2*Fdm));
% determine the frequency spacings of signal after FFT
delta_f = 2*Fdm/N;
% determine time spacing of output samples
delta_T_inv = 1/fs;    %the time duration of a fading waveform.

%   3---Generate complex Gaussian random variables for each of the N/2 positive frequency components of 
%       the noise source.

I_input_time = randn(1,N);
Q_input_time = randn(1,N);
% take FFT
I_input_freq = fft(I_input_time);
Q_input_freq = fft(Q_input_time);
%   4---Construct the negative frequency components of the noise source by conjugating positive frequency 
%       values and assigning these at negative frequency values.

%for i=N/2+1:N
%    I_input_freq(i) = conj( I_input_freq(N+1-i) );
%    Q_input_freq(i) = conj( Q_input_freq(N+1-i) );
%end

%   5---Multiply the in-phase and quadrature noise source by the fading spectrum sqrt(SEZ(f));
SEZ(1) = 1.5/(pi*Fdm);
% 0 < f < Fdm 
for i=2:N/2
   f(i) = (i-1)*delta_f;
   SEZ(i) = 1.5/(pi*Fdm*sqrt(1-(f(i)/Fdm)^2));
   SEZ(N-i+2) = SEZ(i);
end

%SEZ0 = 1.5/(pi*Fm);
% 0 < f < Fm 
%for j=1:(N/2-1)
%   d_f(j) = (j-1/2)*delta_f;
%   SEZ(j) = SEZ0/sqrt(1-(d_f(j)/Fm)^2);
%   SEZ(N-j+1) = SEZ(j);
%end

% use a polynomial fit to get the component at f = Fm(j=N/2 or j=N/2+1)
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 );

I_output_freq = I_input_freq .* sqrt(SEZ);
Q_output_freq = Q_input_freq .* sqrt(SEZ);
%   6---Perform an IFFT to get two N-length time series, 
% take IFFT
I_temp = [I_output_freq(1:N/2) zeros(1,N_inv-N) I_output_freq(N/2+1:N)];
I_output_time = 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 = ifft(Q_temp);

%   7---normalize and compute rms level
rms = sqrt( mean( (I_output_time(1:Ns)+j*Q_output_time(1:Ns)).*conj(I_output_time(1:Ns)+j*Q_output_time(1:Ns))));
imp_rsp_c = (I_output_time(1:Ns)+j*Q_output_time(1:Ns))/rms;
imp_rsp_c=conj(imp_rsp_c');% tran row to column vector
%figure
%plot(1:Ns,abs(imp_rsp_c),'b-');
return;