channel_jake.m

Jakes channel simulation

function [H,t]=channel_jake(DS,Nsamp,M);
%This is a small example illustrates the channel variation with 
%doppler frequncy 
% fd dopper shift in Hz
% Nsamp  = number of slots * length of slot(156)
% M Number of channels needed, ie. number of Tx antennas
dt = 0.001; %Sample time Leave this at 0.001

if nargin == 0, 
M = 2;
Nsamp = 2e3;
DS = 0.01 %Time doppler fading rate
end

L = 2^11; %Size of fft for demo purposes
fd = DS/dt; %Doppler shift (Hz)
%CHOOSE ONE OF THE FOLLOWING !!!
%
%If you want the upper time limit to increase with increasing Nsamp. 
t=dt*(1:Nsamp);

%If you want an increased resolution but use the same upper time limit.
%t=linspace(1e-6,1e-0,Nsamp);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H=rayleigh2(t,fd,M);          % Rayleigh fading
H = conj(H)';


if nargin==0,
    h = H;
    figure(5);clf;
    subplot(2,1,1)
    plot(t,real(h(1,:)));
    tle=['Time-Doppler Shift: ' num2str(fd*dt) '  Sample time ' num2str(dt) '  Doppler (Hz): ' num2str(fd) ];
    title(tle)
    xlabel('Time (sec)')
    ylabel('Amplitude (real component)')
    xlim([0 t(end)])
    subplot(2,1,2);
    plot(imag(h(1,:)));
    xlabel('Number Samples')
    ylabel('Amplitude (imaginary component)')
    xlim([1 Nsamp])
    
%     figure(6); clf;
%     subplot(2,1,1)
%     plot(t,real(h(1,:)));
%     tle=['Time-Doppler Shift: ' num2str(fd*dt) '  Sample time ' num2str(dt) '  Doppler (Hz): ' num2str(fd) ];
%     title(tle)
%     xlabel('Time (sec)')
%     ylabel('Voltage')
%     xlim([0 t(end)])
%     subplot(2,1,2)
%     data = abs(h(1,:));
%     data = data./max(data);
%     semilogy(t,data);
%     xlabel('Time (sec)')
%     ylabel('Power (dB)')
    
    %Fourier analysis signal
    %Real
    h = h(1,:) .* exp(i*2*pi*20*t);
    Yr = fft(real(h),L);
    Pr = Yr.*conj(Yr)/L;
    f = inv(dt)/L*(0:L/2-1);
    %Imaginary
    Yi = fft(imag(h),L);
    Pi = Yi.*conj(Yi)/L;
    f = inv(dt)/L*(0:L/2-1);
    
    figure(7);
    subplot(2,1,1)
    plot(f,Pr(1:L/2))
    title('Power spectral density')
    xlabel('Frequency (Hz)')
    axis([0 5*DS/dt 0 max(Pr(1:L/2))])
    subplot(2,1,2)
    plot(f,Pi(1:L/2))
    xlabel('Frequency (Hz)')
    axis([0 5*DS/dt 0 max(Pi(1:L/2))])


end