S = RandStream('swb2712', 'Seed', 12345); % Set a local random number stream
M = 2; % Modulation order
hModem = modem.pskmod(M); % 2-PSK modulator object
Rsym = 10e3; % Input symbol rate
Rbit = Rsym * log2(M); % Input bit rate
Nos = 4; % Oversampling factor
ts = (1/Rbit) / Nos; % Input sample period
tau = [0 0.4 0.9]*1e-6; % Path delays, in seconds
pdb = [0 -15 -20]; % Average path gains, in dB
dop = doppler.rounded; % Doppler spectrum, with default parameters
fd = 0.5; % Maximum Doppler shift for all paths (identical)
Nt = 2; % Number of transmit antennas
Nr = 1; % Number of receive antennas
rho = 0.7; % Correlation coefficient = antenna correlation
h = mimochan(Nt, Nr, ts, fd, tau, pdb); % MIMO channel object
h.KFactor = 4; % Rician K-factor on first path
h.DopplerSpectrum = dop; % Doppler spectrum of MIMO object
h.TxCorrelationMatrix = toeplitz([1 rho]); % Transmit correlation matrix
% After each frame is processed, the channel is not reset: this is necessary
% to preserve continuity across frames.
h.ResetBeforeFiltering = 0;
% This setting is needed to store the path gains.
h.StorePathGains = 1;
Nsamp = 1.5e6; % Total number of channel samples
Nsamp_f = 1000; % Number of samples per frame
Nframes = Nsamp/Nsamp_f; % Number of frames
out = zeros(Nsamp, Nr);
y1 = zeros(Nsamp, Nt); y2 = zeros(Nsamp, Nt); y3 = zeros(Nsamp, Nt);
for iFrames = 1:Nframes
inputSig = modulate(hModem, randi(S, [0 M-1], Nsamp_f, Nt));
idx = (1:Nsamp_f)+(iFrames-1)*Nsamp_f;
out(idx,:) = filter(h, inputSig);
for it = 1:Nt
% For each transmit antenna, store gains of all three paths
y1(idx,it) = h.PathGains(:,1,it,1);
y2(idx,it) = h.PathGains(:,2,it,1);
y3(idx,it) = h.PathGains(:,3,it,1);
end
end
Hs = spectrum.welch('Hamming', Nsamp/5, 50);
figure;
psd(Hs, y2(:,1), 'Fs', 1/ts, 'SpectrumType', 'twosided', 'Centerdc', true)
axis([-0.1/10 0.1/10 -80 10]);
legend('Simulation');
f = -fd: 0.01 :fd;
a = [1 -1.72 0.785]; % Parameters of the rounded Doppler spectrum
Sd = 1/(2*fd*(a(1)+a(2)/3+a(3)/5)) * ( a(1) + a(2)*(f/fd).^2 + a(3)*(f/fd).^4 );
Sd = Sd * 10^(-15/10); % Scaling by average path power
hold on;
plot(f(Sd>0)/1e3, 10*log10(Sd(Sd>0)), 'k--');
legend('Simulation', 'Theory');
figure;
psd(Hs, y2(:,2), 'Fs', 1/ts, 'SpectrumType', 'twosided', 'Centerdc', true)
axis([-0.1/10 0.1/10 -80 10]);
legend('Simulation');
hold on;
plot(f(Sd>0)/1e3, 10*log10(Sd(Sd>0)), 'k--');
legend('Simulation', 'Theory');
figure;
semilogy(abs(y1(:,1)),'b');
hold on; grid on;
semilogy(abs(y1(:,2)),'r');
legend('First transmit link', 'Second transmit link');
title('Fading envelopes for two transmit links of Path 1');
figure;
semilogy(abs(y2(:,1)),'b');
hold on; grid on;
semilogy(abs(y2(:,2)),'r');
legend('First transmit link', 'Second transmit link');
title('Fading envelopes for two transmit links of Path 2');
figure;
semilogy(abs(y3(:,1)),'b');
hold on; grid on;
semilogy(abs(y3(:,2)),'r');
legend('First transmit link', 'Second transmit link');
title('Fading envelopes for two transmit links of Path 3');
TxCorrMatrixPath1 = corrcoef(y1(:,1),y1(:,2)).'
TxCorrMatrixPath2 = corrcoef(y2(:,1),y2(:,2)).'
TxCorrMatrixPath3 = corrcoef(y3(:,1),y3(:,2)).'