% Code to demonstrate the implementation of a narrowband fading channel
% with memory by use of a correlation matrix for time based filtering
% detailed in section 4.1.2
% Tim Brown, January 2012

N = 1000; % Arbritrary number of samples, 1000 is useful to see clear comparison.
Delta_t = 0.1; % Arbritrary change in time in seconds

% Create the rows of the correlation matrix, in this case defining it based
% on the Bessel function suited to Rayleigh channels
Phase_range = 2.*pi.*[0:(60/(N-1)):60]; % Representative of moving 30 wavelengths
Delta_phase = 2*pi*10/(N-1); % Change in phase per sample
R = zeros(N,N); % Forming the full correlation matrix as in the book (though not fully used)
for nn = 1:N 
    R(nn,:) = besselj(0,Phase_range-((nn-1)*Delta_phase));
end
% Create random Rayleigh samples without memory
h = (randn(N, 1) + j*randn(N, 1))/ sqrt(2);
% Perform correlation on the samples to create channel with memory
row1= R(1,:); % Note this line is different to as in the book
column1 = row1';
matrixR = toeplitz(row1, column1);
W = sqrtm(matrixR);
w=W*h;

% Take the Fourier transform of both the channel with and without memory
fft_h = fftshift(fft(h));
fft_w = fftshift(fft(w));

% Plot a comparison of the channels with and without memory in the
% frequency domain
f = [-(0.5/Delta_t):((1/Delta_t)/(N-1)):(0.5/Delta_t)];
figure
plot(f,20.*log10(abs(fft_h)),f,20.*log10(abs(fft_w)))
set(gca,'FontSize',10)
set(gcf,'color','white')
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
legend('With Memory','Without Memory')

% Plot in the time domain, the first 100 for clarity, though this can be
% changed by adjusting 'range' below
range = 1:100;
t = (0:(N-1))*Delta_t;
figure
plot(t(range),20.*log10(abs(h(range))),t(range),20.*log10(abs(w(range))))
set(gca,'FontSize',10)
set(gcf,'color','white')
xlabel('Time (s)')
ylabel('Magnitude (dB)')
legend('With Memory','Without Memory')