rayleigh1.m

simulation of rayleigh channel

function [c]=rayleigh(N,v,fc,fs)

% I.i.d. Gaussian samples from which fading process will be generated
N= 100000; %numberOfChannelSamples
iidGaussian=(1/sqrt(2))*(randn(1,100000)+sqrt(-1)*randn(1,100000)); 
%hist(iidGaussian)

%CREATE UNIT POWER COMPLEX IID GAUSSIAN

% Compute the doppler frequency as a function of mobile velocity and
% % carrier frequency
% v = 30; % [km/hr]
v=  30*5/18; % [m/s]
 fc = 2e9 ; % PCS band
% lambda = (3*10^8)/2e9; %wavelength
 doppler = fc*v/(3*10^8)  % maximum doppler frequency [Hz]
% 
fprintf('Doppler frequency: %.2f Hz\n', doppler);

% Fading process oversampling rate
samplingFreq = 2 * 2 * doppler       % 2x oversampled
samplingPeriod =  1/samplingFreq; % Sampling Period  

% Fading process will be created by filtering i.i.d. Gaussian 
% with a filter with Clark's spectrum

% Compute the autocorrelation function of fading process:
% J0(2 pi f_d tau)
filterLength = 2^10;
timeRange = ( (1:filterLength)-filterLength/2 ) * samplingPeriod;


% Create the spectrum of fading process
FilterSpectrum = [];
freqRange = linspace(-samplingFreq/2, samplingFreq/2, filterLength);
for freq = freqRange;
    if (abs(freq) < doppler)
        FilterSpectrum = [FilterSpectrum (1.5)/(pi*doppler*sqrt(1-(freq/doppler)^2))];
    else
        FilterSpectrum = [FilterSpectrum 0];
    end
end
figure(1)
plot(freqRange/doppler,FilterSpectrum), xlabel('Relative Frequency f/fd'),ylabel('Magnitude');
title('clarks model Power spectrum');

% Do a cyclic shift on filter Spectrum so that the zero frequency approximately coincides with the
% first entry of the array. 
FilterSpectrum = [FilterSpectrum(filterLength/2+1:end) fliplr(FilterSpectrum(filterLength/2+1:end))];

% Compute the time-domain filter with which i.i.d. Gaussian will be
% filtered- use ifft function
FilterFreq = sqrt(FilterSpectrum);
FilterTimeDomain =ifft(FilterFreq);
FilterTimeDomain = FilterTimeDomain/norm(FilterTimeDomain);
FilterTimeDomain = real(FilterTimeDomain); % Eliminate spurious complex entries that may arise 
rayleighProcess=fftfilt(FilterTimeDomain,[iidGaussian zeros(1,filterLength)]);
rayleighProcess=rayleighProcess(1+filterLength/2:end-filterLength/2);            
[xryf,yryf]=fCDF(abs(rayleighProcess)); 
figure(2)
plot(xryf,yryf),xlabel('Magnitude of complex envelope');
ylabel('Probability the abscissa is not exceeded'),title('CDF of Rayleigh series after filtering');
%resample the rayleigh process to match the input signal sampling rate
[p,q]=rat(fs/samplingFreq,0.0001);
c=resample(rayleighProcess,p,q);
% % Compute the Rayleigh process by convolving i.i.d. Gaussian with the
% % computed filter- use conv function
% hist(abs(rayleighProcess))
% sig=rand(1,4024)>0.5;
% sig=reshape(sig,length(sig)/2,2);
% sig=bi2de(sig,'left-msb')';
% j=sqrt(-1);
% sig_mod=[];
% for i=1:2012
%     if sig[i]==1
%         gsig_mod(i)= -1+j;
%     elseif sig[i]==0
%         sig_mod(i)=-1-j;
%     elseif sig[i]==2
%         sig_mod(i)=1-j;
%     else
%         sig_mod(i)=1+j;
%     end
% end