Rayleigh fading simulation

[irpotential]

Hi, I've search around a lot and cannot figure out what I'm doing wrong here.

The problem states to generate a time sequence of 8192 sample values of a Rayleigh fading signal for a doppler frequency of 20 Hz.

The process we have been told to use is outlined as follows:
1) The PSD of our received signal r(t) with average power P_r is S_r=P_r/4pi*f_D * 1/√1-(f/fD)²
2) Determine sample time interval, spacing between frequency components, etc.
3) Generate complex Gaussian random variables for each of the N/2 positive frequency components of the noise source: Use different sets of random variables for the in-phase and quadrature components
4) Construct the negative components of the noise source by conjugating positive frequency values and then assigning them to the corresponding negative frequency components
5) Multiply the in-phase and quadrature noise sources by the fading spectrum √S_r
6) Perform IFFT on the resulting frequency domain signals from the in-phase and quadrature arms to get two time series of length N. Then add the squares of each signal point in time to create an N-point time series. Note that each quadrature arm should be a real signal after taking the IFFT.
7) Take the square root of the sum from 6 to obtain an N-point time series r(t), which simulates the Rayleigh fading signal with the proper Doppler spread and time correlation

MATLAB Code:


clear all; close all; clc;

% data given in the problem
f_D = 20; % doppler freq (Hz)
N = 8192; % number of points

% data not given but needed
P_r = 1; % power received (W)

% sample frequency
f_s = 2*f_D/N;

% frequency for each point
% NOTE: S_r(f_D) = inf
f = f_s:f_s:f_D-f_s;

% PSD of r(t)
S_r = P_r./(4*pi*f_D*sqrt(1-(f/f_D).^2));

% create the gaussian random variables
A = complex(randn(1, N/2-1),randn(1, N/2-1));
B = complex(randn(1, N/2-1),randn(1, N/2-1));

% multiply the gaussian variables by the fading spectrum, sqrt(S_r)
A = A.*sqrt(S_r);
B = B.*sqrt(S_r);

%now make the negative frequency components
A = [fliplr(conj(A)) A];
B = [fliplr(conj(B)) B];

% take the ifft and square each value
a = ifft(A).^2;
b = ifft(B).^2;

% compute r
r = sqrt(a+b);

plot(20*log10(r))


please help!