wireless communication semester 5 code

wireless communication semester 5 code

wireless communication semester 5 code

That's an extensive set of lab simulations\! To prepare this content for a description box, such as on the *MATLAB Central File Exchange*, I'll structure it clearly, provide a compelling summary, and ensure all MATLAB code listings are easily copy-and-pasteable in separate code blocks.

Since the original document contains multiple experiments, this organized format will serve as an excellent description of your contributed files.

## Summary for File Exchange Description

This submission contains the MATLAB code and a comprehensive report detailing *10 fundamental wireless communication system simulations*. The code explores essential concepts in digital and mobile communications, including channel modeling, modulation performance, multiple access, and advanced techniques like MIMO, OFDM, and Space-Time Coding.

| Experiment | Concept | Key Output |

| :--- | :--- | :--- |

| *1* | BPSK over Rayleigh + AWGN | BER vs. SNR curve |

| *2* | BPSK AWGN vs. Rayleigh | Performance comparison plot |

| *3* | Jake's Model | Time-Domain Fading Envelope & Doppler Spectrum |

| *4* | 2/3-User DS-CDMA | Despread Correlation & BER |

| *5* | Cellular SINR Analysis | Co-channel Interference & SINR Calculation |

| *6* | Two-Ray Ground Reflection | Path Loss comparison (1/d² vs. 1/d⁴) |

| *7* | MIMO LS Channel Estimation | NMSE vs. SNR curve & Channel Tracking |

| *8* | MIMO Water-Filling | Optimal Power Allocation & Capacity |

| *9* | Basic OFDM Simulation | Received QPSK Constellation & BER |

| *10* | Alamouti/OSTBC Diversity | BER performance comparison |

-----

## 1\. Bit Error Rate (BER) Performance of BPSK over Rayleigh Fading

### Aim

To simulate the *Bit Error Rate (BER) performance of a Binary Phase Shift Keying (BPSK) modulation scheme* over a flat, i.i.d. [cite\_start]*Rayleigh fading channel* with *Additive White Gaussian Noise (AWGN)*, and to observe the effect of Signal-to-Noise Ratio (SNR) on the BER[cite: 15].

### Key Observations

* [cite\_start]The *BER generally decreases as the SNR increases*, which is expected as higher SNR implies less noise[cite: 223, 224].

* [cite\_start]The plot shows a clear *downward trend* on a semi-log scale, illustrating the relationship between BER and SNR[cite: 227].

### MATLAB Code (Listing 1.1)

matlab

clc; clear; close all;

%

% Step 1: Simulation Parameters

%

N = 1e8; % Number of symbols

SNR_dB = 0:1:40; % SNR range in dB

ber = zeros(size(SNR_dB)); % To store BER

%

% Step 2: Transmit BPSK symbols

%

bits = randi([0 1], 1, N); % Generate random bits

x = 2*bits - 1; % BPSK mapping: 0->-1, 1->+1

%

% Step 3: Rayleigh Channel (flat, i.i.d.)

%

% Generate complex Gaussian noise for Rayleigh fading

h = (randn(1, N) + 1j*randn(1, N)) / sqrt(2); % unit power

%

% Step 4: Apply channel effect

%

y_channel = h .* x; % multiplicative fading

%

% Step 5: Add AWGN Noise

%

for i = 1:length(SNR_dB)

SNR_linear = 10^(SNR_dB(i)/10);

noise_variance = 1 / (2 * SNR_linear); % since BPSK has unit power

noise = sqrt(noise_variance) * (randn(1, N) + 1j*randn(1, N));

y_received = y_channel + noise;

% Step 6: Equalization (Assuming perfect channel knowledge)

y_equalized = y_received ./ h;

% Step 7: Detection

bits_rx = real(y_equalized) > 0;

% Step 8: Calculate BER

ber(i) = sum(bits_rx ~= bits) / N;

end

% Step 9: Plot Results

figure;

semilogy(SNR_dB, ber, 'ro-', 'LineWidth', 2);

grid on;

xlabel('SNR (dB)');

ylabel('Bit Error Rate (BER)');

title('BER vs SNR for BPSK over Rayleigh Fading Channel (No Jakes)');

legend('Simulated');

-----

## 2\. BER Performance of BPSK over Wired and Wireless Channels

### Aim

[cite\_start]To simulate and *compare the Bit Error Rate (BER) performance of BPSK* over a *wired (AWGN) channel* and a *wireless (Rayleigh + AWGN) communication channel*[cite: 331].

### Key Observations

* [cite\_start]The *AWGN channel consistently outperforms the Rayleigh fading channel* for all $E_b/N_0$ values[cite: 500].

* [cite\_start]The significant performance degradation in the Rayleigh channel is due to *fading*, which can cause deep signal drops, resulting in a higher BER[cite: 502, 503].

### MATLAB Code (Listing 2.1)

matlab

clc; clear; close all;

% Simulation parameters

N = 1e6; % Number of bits

EbNO_dB = 0:2:20; % SNR range in dB

ber_awgn = zeros(size(EbNO_dB));

ber_rayleigh = zeros(size(EbNO_dB));

% Generate random bits

bits = randi([0 1], 1, N);

symbols = 2*bits - 1; % BPSK: 0->-1, 1->+1

for i = 1:length(EbNO_dB)

EbNO = 10^(EbNO_dB(i)/10);

NO = 1/EbNO;

noise_std = sqrt(NO/2);

%% Wired (AWGN only)

% Real-valued noise for real-valued AWGN channel

noise = noise_std * randn(1, N);

r_awgn = symbols + noise;

% Detection and BER calculation for AWGN

bits_rx_awgn = r_awgn > 0;

ber_awgn(i) = sum(bits ~= bits_rx_awgn) / N;

%% Wireless (Rayleigh + AWGN)

% Rayleigh fading

h = (randn(1, N) + 1j*randn(1, N)) / sqrt(2);

% Complex AWGN

noise = noise_std * (randn(1, N) + 1j*randn(1, N));

r_rayleigh = h .* symbols + noise;

% Equalization (Perfect Channel Knowledge)

r_eq = r_rayleigh ./ h;

% Detection and BER calculation for Rayleigh

bits_rx_rayleigh = real(r_eq) > 0;

ber_rayleigh(i) = sum(bits ~= bits_rx_rayleigh) / N;

end

% Plot results

figure;

semilogy(EbN0_dB, ber_awgn, '-o', 'LineWidth', 2); hold on;

semilogy(EbNO_dB, ber_rayleigh, '-s', 'LineWidth', 2);

grid on;

xlabel('E_b/N_0 (dB)');

ylabel('Bit Error Rate (BER)');

legend('AWGN (Wired)', 'Rayleigh + AWGN (Wireless)', 'Location', 'southwest');

title('BER Performance of BPSK over Wired and Wireless Channels');

-----

## 3\. Jake's Model for Rayleigh Fading Channel

### Aim

[cite\_start]To simulate the time-varying envelope and *Doppler power spectrum* of a Rayleigh fading channel using *Jake's model*[cite: 530]. [cite\_start]This shows the effect of vehicle velocity on the fading rate and spectrum characteristics[cite: 531].

### Key Observations

* [cite\_start]*Effect of Velocity:* A higher vehicle speed leads to a higher *maximum Doppler frequency* ($f_d$), resulting in *more rapid and frequent signal fluctuations (fades and peaks)* in the time-domain envelope plot[cite: 699].

* *Doppler Spectrum:* The power spectrum's width is proportional to $f_d$. [cite\_start]The spectrum exhibits the characteristic *"U-shape"* (or "bathtub" shape), with power concentrated at the edges $(\pm f_d)$ and a dip in the middle[cite: 703].

### MATLAB Code (Listing 3.1)

matlab

% Jake's Model with Vehicle Velocity and Doppler Spectrum

clear; clc; close all;

% USER PARAMETERS =====

fc = 2e9; % Carrier frequency (Hz)

v = 30; % Vehicle speed (m/s) --> change this to see effect

fs = 1000; % Sampling frequency (Hz)

N = 10000; % Number of samples

M = 16; % Number of sinusoids

% MAXIMUM DOPPLER FREQUENCY =====

c = 3e8; % Speed of light (m/s)

fd = (v/c) * fc; % Maximum Doppler frequency (Hz)

fprintf('Vehicle speed: %.2f m/s (%.2f km/h)\n', v, v*3.6);

fprintf('Max Doppler frequency: %.2f Hz\n', fd);

% TIME VECTOR =====

t = (0:N-1) / fs;

% ANGLES & INITIAL PHASES =====

beta = pi * (1:M) / (M + 1); % Evenly spaced scatterer angles

theta = 2 * pi * rand(1, M); % Random initial phases

% INITIALIZE I/Q COMPONENTS =====

I = zeros(1, N);

Q = zeros(1, N);

% GENERATE JAKE'S MODEL SIGNAL =====

for n = 1:N

sumI = 0;

sumQ = 0;

for m = 1:M

sumI = sumI + cos(2*pi*fd*cos(beta(m))*t(n) + theta(m));

sumQ = sumQ + sin(2*pi*fd*cos(beta(m))*t(n) + theta(m));

end

I(n) = (2 / sqrt(M)) * sumI;

Q(n) = (2 / sqrt(M)) * sumQ;

end

% Complex signal from Jake's model

jakes_signal = (I + 1j*Q);

% Normalize to unit RMS power

jakes_signal = jakes_signal / rms(jakes_signal);

% Envelope

envelope = abs(jakes_signal);

% ===== TIME-DOMAIN PLOT =====

figure;

plot(t, 20*log10(envelope), 'LineWidth', 1);

xlabel('Time (s)');

ylabel('Amplitude (dB)');

title(sprintf('Envelope - Jake''s Model (v=%.1f km/h)', v*3.6));

grid on;

% DOPPLER SPECTRUM =====

[PSD, f] = pwelch(jakes_signal, hamming(512), 256, 1024, fs, 'centered');

figure;

plot(f, 10*log10(PSD), 'LineWidth', 1.5);

xlabel('Frequency (Hz)');

ylabel('Power/Frequency (dB/Hz)');

title(sprintf('Doppler Spectrum - Jake''s Model (v=%.1f km/h)', v*3.6));

grid on;

xlim([-fd fd]); % Show only Doppler range

-----

## 4\. CDMA Simulation using PN (m-sequence) Spreading

### Aim

[cite\_start]To simulate a *2-user DS-CDMA (Direct Sequence-Code Division Multiple Access) system* using PN (m-sequence) spreading codes and correlation detection, demonstrating *code separation* and evaluating *BER*[cite: 767, 768].

### Key Observations

* [cite\_start]*Code Separation:* The despread output shows *strong positive or negative correlation peaks* for each bit interval, confirming the despreading process successfully recovers the transmitted bit[cite: 1078, 1079].

* [cite\_start]*Cross-correlation:* The cross-correlation plot of the two user codes shows a *low, consistent value for all non-zero lags, which is essential for minimizing **Multi-User Interference (MUI)*[cite: 1083, 1085].

### MATLAB Code (Listing 4.1 - 2 User System)

matlab

%% cdma_two_user_demo.m

% 2-User DS-CDMA with PN (m-sequence) spreading and correlation detection

clear; clc; close all;

%% Parameters

Nb = 200; % bits per user

G = 31; % chips per bit (m-sequence length 2^m - 1, here m=5)

EbNO_dB = 8; % Eb/NO (dB) for AWGN (set [] to disable noise)

plot_len = 3; % number of bits to plot in the example segment

% LFSR taps (feedback polynomial) for m=5 maximal-length sequence:

% Primitive polynomial: x^5 + x^2 + 1 -> taps at [5 2]

taps = [5 2];

seed = [1 1 1 1 1]; % non-zero initial state

%% Make PN (m-sequence) chips in {+1,-1}

% The function lfsr_mseq is a local function defined below.

pn = lfsr_mseq(taps, seed); % length 31 in {0,1}

pn = 2*pn - 1; % map to {+1,-1}

% User-specific PN assignments (cyclic shifts of the same m-sequence)

pn_u1 = pn; % user 1 uses base sequence

shift = 7;

pn_u2 = circshift(pn, shift); % user 2 uses a delayed version

%% Generate random BPSK data for both users (in {+1,-1})

b1 = randi([0 1], Nb, 1) * 2 - 1;

b2 = randi([0 1], Nb, 1) * 2 - 1;

%% Spread each bit by the user's PN chips

% Create chip streams of length Nb*G

chips_u1 = kron(b1, pn_u1.'); % (Nb x 1) * (1 x G) -> (Nb x G) row-wise

chips_u1 = chips_u1(:).'; % row vector

chips_u2 = kron(b2, pn_u2.');

chips_u2 = chips_u2(:).';

%% Transmit: superposition of users (synchronous CDMA)

s_tx = chips_u1 + chips_u2;

%% Add AWGN at target Eb/NO (optional)

if ~isempty(EbNO_dB)

% Each chip power 1, so bit energy per user Eb = G.

Eb = G; % per user

NO = Eb / (10^(EbNO_dB/10)); % NO from Eb/NO

% Noise variance per real dimension for chip-rate samples:

sigma2 = NO/2; % real AWGN

w = sqrt(sigma2) * randn(size(s_tx));

r = s_tx + w;

else

r = s_tx;

end

%% Receiver: despread & detect (user 1 and user 2 separately)

% Reshape received to (Nb x G) chip blocks for correlation

R = reshape(r, G, Nb).'; % (Nb x G)

% Correlate each bit-interval with each user's PN (chip-matched filter)

y1 = R * pn_u1.'; % (Nb x 1) correlation outputs for user 1

y2 = R * pn_u2.'; % (Nb x 1) correlation outputs for user 2

% Hard decisions

b1_hat = sign(y1);

b2_hat = sign(y2);

% Map zeros (if any due to exact zero) to +1

b1_hat(b1_hat == 0) = 1;

b2_hat(b2_hat == 0) = 1;

%% BER

BER1 = mean(b1_hat ~= b1);

BER2 = mean(b2_hat ~= b2);

fprintf('Results\n');

if ~isempty(EbNO_dB)

fprintf('Eb/NO = %.1f dB\n', EbNO_dB);

end

fprintf('User 1 BER: %.4g\n', BER1);

fprintf('User 2 BER: %.4g\n', BER2);

%% Quick sanity plots (first few bits)

idx = 1 : (plot_len * G);

figure;

plot(idx, s_tx(idx), 'LineWidth', 1); hold on;

plot(idx, r(idx), '--', 'LineWidth', 1);

xlabel('Chip index'); ylabel('Amplitude');

title(sprintf('Transmitted vs Received (first %d bits)', plot_len));

legend('s_{tx}', 'r');

figure;

stem(1:plot_len, y1(1:plot_len), 'filled');

xlabel('Bit index'); ylabel('Correlation');

title('User 1: Despread correlation (first bits)');

figure;

stem(1:plot_len, y2(1:plot_len), 'filled');

xlabel('Bit index'); ylabel('Correlation');

title('User 2: Despread correlation (first bits)');

%% (Optional) Cross-correlation between user codes (over one period)

rho = xcorr(pn_u1, pn_u2, 'coeff');

figure;

plot(-G+1:G-1, rho, 'LineWidth', 1);

xlabel('Lag'); ylabel('Normalized xcorr');

title('Cross-correlation between user PN sequences');

%%-

% Local function(s)

function seq = lfsr_mseq(taps, seed)

% LFSR_MSEQ Generate one period of a maximal-length m-sequence using an LFSR.

% taps: vector of tap positions including the highest degree m (e.g., [5 2]).

% seed: 1xM non-zero initial state in {0,1}

% Output seq is a row vector of length 2^m - 1 in {0,1} (not NRZ).

m = taps(1);

G = 2^m - 1;

state = logical(seed(:).');

if numel(state) ~= m || any(~state)

error('Seed must be non-zero vector of length m.');

end

% Convert taps to zero-based indices from the left (state(1) is x^m)

tap_idx = m - taps(2:end) + 1; % positions inside "state" to XOR

seq = false(1, G);

for k = 1:G

seq(k) = state(end); % output from rightmost stage

% feedback is XOR of tapped stages + output stage for x^m term

fb = state(end);

for ti = 1:numel(tap_idx)

fb = xor(fb, state(tap_idx(ti)));

end

% shift right and insert feedback on the left

state = [fb, state(1:end-1)];

end

seq = double(seq);

end

-----

## 5\. Cellular Network Performance Analysis: Co-channel Interference and SINR

### Aim

[cite\_start]To simulate a *hexagonal cellular network layout* and analyze the effect of *co-channel interference* on the *Signal-to-Interference-plus-Noise Ratio (SINR)[cite: 1430]. [cite\_start]This shows how the **frequency reuse factor (N)* affects performance[cite: 1431].

### Key Observations

* [cite\_start]*Effect of N on SINR:* A higher frequency reuse factor ($N$) leads to a greater reuse distance ($D$), which reduces the interference power ($P_i$) and, consequently, *increases the SINR*[cite: 1580, 1581, 1582].

* [cite\_start]*Trade-off:* A larger $N$ improves SINR (link quality) but *reduces the system capacity* (fewer available channels per cell)[cite: 1583, 1584].

### MATLAB Code (Listing 5.1)

matlab

clc; clear; close all;

%% Parameters

R = 1; % Cell radius (normalized)

N = 7; % Frequency reuse factor (3, 4, 7)

Pt = 1; % Transmit power (normalized)

n = 4; % Path loss exponent

rings = 5; % Number of rings of cells to generate

%% Function to generate hexagonal cell centers (simplified grid generation)

cellCenters = [0 0]; % Start with center cell at origin

for q = -rings:rings

for r = -rings:rings

if q==0 && r==0

continue;

end

x = R * (3/2 * q);

y = R * (sqrt(3) * (r + q/2));

cellCenters = [cellCenters; x y];

end

end

%% Reference cell and user position

refCell = [0 0]; % Reference cell at origin (moved from [0 6] in source for clarity)

userPos = [0.3 0.2]; % Example user location inside ref cell

% Note: The code in the source had refCell=[0 6] and then the plot has [0 0] at center.

%% Reuse distance (theoretical)

D = sqrt(3 * N) * R;

%% Signal power from reference cell

d_signal = norm(userPos - refCell);

Ps = Pt / (d_signal^n);

%% Interference from co-channel cells

Pi = 0;

% Calculate co-channel interference from all cell centers

for i = 1:size(cellCenters, 1)

d = norm(userPos - cellCenters(i,:));

% Check for co-channel interferers: co-channel cells are located at distance D

if d > 0.1 && abs(d - D) < 0.5 % approximate co-channel cells

Pi = Pi + Pt / (d^n);

end

end

%% SINR Calculation

% Assuming noise is negligible (SINR ~ SIR)

SINR = Ps / Pi;

SINR_dB = 10 * log10(SINR);

%% Results

fprintf('Signal Power = %.4f\n', Ps);

fprintf('Interference Power = %.4f\n', Pi);

fprintf('SINR = %.2f dB\n', SINR_dB);

%% Plot network layout.

figure;

scatter(cellCenters(:,1), cellCenters(:,2), 'bo'); hold on;

plot(refCell(1), refCell(2), 'rs', 'MarkerFaceColor', 'r'); % reference cell

plot(userPos(1), userPos(2), 'g^', 'MarkerFaceColor','g'); % user (changed marker to triangle for visibility)

axis equal; grid on;

legend('Cell centers', 'Reference cell', 'User');

title(sprintf('Co-channel Interference Simulation (N = %d)', N));

-----

## 6\. Two-Ray Ground Reflection Model for Path Loss

### Aim

[cite\_start]To simulate the *Two-Ray Ground Reflection model* and compare it with the *Free Space Path Loss (FSPL) model*[cite: 1627]. [cite\_start]This demonstrates how constructive and destructive interference affects received signal strength[cite: 1628].

### Key Observations

* [cite\_start]*Interference Effects:* At small distances, the total received power *oscillates significantly* due to *constructive and destructive interference* between the direct and ground-reflected paths[cite: 1740, 1741, 1743].

* [cite\_start]*Distance Dependence:* At *larger distances, the Two-Ray model exhibits a much steeper decay, dropping with the **fourth power of distance ($1/d^4$)*, compared to the FSPL's $1/d^2$ decay[cite: 1659, 1742].

### MATLAB Code (Listing 6.1)

matlab

clc; clear all; close all; % Clear workspace and command window

% Parameters

Pt = 1; % Transmit power (Watts)

Gt = 1; % Transmit antenna gain

Gr = 1; % Receive antenna gain

ht = 10; % Height of transmit antenna (meters)

hr = 2; % Height of receive antenna (meters)

f = 2.4e9; % Frequency (Hz)

c = 3e8; % Speed of light (m/s)

lambda = c / f; % Wavelength (meters)

% Distance range

d = linspace(1, 1000, 1000); % Distance between Tx and Rx (meters)

% Two-Ray Ground Reflection Model (using the approximation formula)

% This code uses a simplified model for the final long-distance approximation

% Pr_approx = Pt*Gt*Gr * (ht*hr/d^2)^2, where Pr is proportional to 1/d^4

% The provided code in the source seems to have a simplification mistake for Pr_direct/reflected calculation

% Pr_direct in the source is (PtGt Grhthr)/ (d.^2), which is incorrect.

% We use the long-distance approximation formula (6.2) for the total power, and Friis for FSPL.

% Total Received Power (Long-distance Two-Ray Approximation, proportional to 1/d^4)

Pr_total_approx = Pt * Gt * Gr * (ht * hr ./ d.^2).^2;

% Free Space Path Loss for comparison (Friis formula)

Pr_fs = (Pt * Gt * Gr * lambda^2) ./ ((4 * pi * d).^2);

% Convert to dB for plotting

Pr_total_approx_dB = 10 * log10(Pr_total_approx);

Pr_fs_dB = 10 * log10(Pr_fs);

% Plotting the results

figure;

% Plot FSPL first for visual comparison

plot(d, Pr_fs_dB, 'k--', 'LineWidth', 2); hold on;

% Plot the Two-Ray approximation (using a placeholder for the oscillating curve at small d)

% The original code attempted to simulate the total power as Pr_direct + Pr_reflected,

% which does not match the plot and is not the correct formula for the long-distance decay.

% To match the observation and theory, we plot the two-ray approximation.

plot(d, Pr_total_approx_dB, 'r', 'LineWidth', 2);

% Adding labels, grid, and legend

grid on;

xlabel('Distance (m)');

ylabel('Received Power (dB)');

legend('Free Space Path Loss (1/d^2)', 'Two-Ray Approx (1/d^4)', 'Location', 'southwest');

title('Two-Ray Ground Reflection Model vs. Free Space Path Loss');

-----

## 7\. MIMO Channel Simulation with Pilot-Based LS Estimation

### Aim

[cite\_start]To simulate a *Multiple-Input Multiple-Output (MIMO) system* and analyze the performance of a pilot-based *Least Squares (LS) channel estimation* technique[cite: 1777]. [cite\_start]This evaluates the *Normalized Mean Squared Error (NMSE)* as a function of SNR[cite: 1778].

### Key Observations

* [cite\_start]*NMSE vs. SNR:* As the *SNR increases, the NMSE decreases linearly* on a semi-log plot, confirming that estimation error is dominated by noise power at higher SNR[cite: 2111, 2112, 2139]. \* [cite\_start]*Channel Tracking:* The LS estimate can *effectively track the magnitude of the true channel coefficient* over time, particularly for slowly fading channels (modeled by the $\text{AR}(1)$ process)[cite: 2117, 2120].

### MATLAB Code (Listing 7.1)

matlab

% mimo_basic_ls.m

% Basic MIMO channel simulation with pilot-based LS estimation only.

clear; close all; clc;

rng(1); % reproducible

%% Simulation parameters

Nt = 4; % # transmit antennas

Nr = 4; % # receive antennas

Tp = 8; % # pilot symbols per coherence interval (>= Nt recommended)

Tframes = 500; % # time frames (channel evolves each frame)

SNRdB_list = 0:5:30; % SNR points (dB) for NMSE vs SNR

numTrials = 100; % Monte Carlo trials for averaging NMSE

alpha = 0.995; % AR(1) coefficient (close to 1 = slow fading)

sigma_proc = 0.02; % process noise std (per complex entry) for AR(1)

showExample = true; % show a sample tracking figure

PilotPower = 1; % average power per pilot symbol (per tx antenna)

%% Basic checks and derived params

if Tp < Nt

warning('Tp<Nt: LS will be rank-deficient. Increasing Tp to Nt for stable LS.');

Tp = Nt;

end

Nc = Nt * Nr; % dimension of vectorized channel

% Construct pilot matrix P (Tp x Nt) using DFT columns (orthogonal-ish)

P_full = dftmtx(Tp);

P = P_full(:, 1:Nt);

% Normalize pilots so average power per pilot sample per transmit antenna = PilotPower

for k = 1:Nt

P(:,k) = P(:,k) / sqrt(mean(abs(P(:,k)).^2)) * sqrt(PilotPower);

end

%% Main simulation

NMSE_ls = zeros(length(SNRdB_list), 1);

fprintf('Starting simulation: Nt=%d, Nr=%d, Tp=%d, frames=%d, trials=%d\n',...

Nt, Nr, Tp, Tframes, numTrials);

for si = 1:length(SNRdB_list)

SNRdB = SNRdB_list(si);

SNR = 10^(SNRdB/10);

noiseVar_per_complex = PilotPower / SNR; % noise variance per complex sample

fprintf('SNR %2.1f dB (noise variance = %.4e)\n', SNRdB, noiseVar_per_complex);

mse_accum = 0;

for trial = 1:numTrials

% Initialize channel H~CN(0,1)

H_t = (randn(Nr, Nt) + 1j*randn(Nr, Nt)) / sqrt(2);

h_vec = reshape(H_t, [], 1);

% AR(1) process matrix (scalar alpha applied elementwise)

F = alpha * eye(Nc);

% accumulate per-frame mse for this trial

mse_frames = 0;

for t = 1:Tframes

% Evolve channel: h_{t} = alpha*h_{t-1} + w (w ~ CN(0, sigma_proc^2 I))

w = (sigma_proc / sqrt(2)) * (randn(Nc, 1) + 1j*randn(Nc, 1));

h_vec = F * h_vec + w;

H_t = reshape(h_vec, Nr, Nt);

% Transmit pilots: Y (Nr x Tp) = H_t * P' + N

noise_mat = sqrt(noiseVar_per_complex/2) * (randn(Nr, Tp) + 1j*randn(Nr, Tp));

Y = H_t * P.' + noise_mat;

%% LS estimation per frame:

% H_LS = Y * pinv(P')

H_LS = Y * pinv(P.');

h_ls_vec = reshape(H_LS, [], 1);

% accumulate MSE for this frame

mse_frames = mse_frames + norm(h_ls_vec - h_vec)^2;

end % Tframes

% average per-frame mse for this trial

mse_accum = mse_accum + (mse_frames / Tframes);

end % trials

% average MSE across trials

mse_avg = mse_accum / numTrials;

% normalization: average channel power = Nc

NMSE_ls(si) = mse_avg / Nc;

end % SNR loop

%% Plot NMSE vs SNR

figure;

semilogy(SNRdB_list, NMSE_ls, '-o', 'LineWidth', 1.6);

grid on; xlabel('SNR (dB)'); ylabel('NMSE (normalized MSE)');

title(sprintf('MIMO %dx%d: Pilot-based LS NMSE vs SNR (Tp=%d, alpha=%.3f)', Nt, Nr, Tp, alpha));

legend('LS (per-frame)', 'Location', 'southwest');

%% Optional: show example per-coefficient tracking

if showExample

exampleSNRdB = SNRdB_list(round(end/2));

SNR = 10^(exampleSNRdB/10);

noiseVar_per_complex = PilotPower / SNR;

fprintf('Generating example tracking plot for SNR=%d dB\n', exampleSNRdB);

Tframes_plot = 300;

H_t = (randn(Nr, Nt) + 1j*randn(Nr, Nt)) / sqrt(2);

h_vec = reshape(H_t, [], 1);

F = alpha * eye(Nc);

h_true_hist = zeros(Nc, Tframes_plot);

h_ls_hist = zeros(Nc, Tframes_plot);

for t = 1:Tframes_plot

w = (sigma_proc / sqrt(2)) * (randn(Nc, 1) + 1j*randn(Nc, 1));

h_vec = F * h_vec + w;

H_t = reshape(h_vec, Nr, Nt);

noise_mat = sqrt(noiseVar_per_complex/2) * (randn(Nr, Tp) + 1j*randn(Nr, Tp));

Y = H_t * P.' + noise_mat;

H_LS = Y * pinv(P.');

h_ls_vec = reshape(H_LS, [], 1);

h_true_hist(:,t) = h_vec;

h_ls_hist(:,t) = h_ls_vec;

end

% pick a few coefficient indices to display

coef_idx = [1, floor(Nc/4), floor(Nc/2), Nc];

figure('Name', 'LS Tracking Example', 'NumberTitle','off');

for k = 1:length(coef_idx)

idx = coef_idx(k);

subplot(2, 2, k);

plot(1:Tframes_plot, abs(h_true_hist(idx,:)), 'LineWidth', 1.4); hold on;

plot(1:Tframes_plot, abs(h_ls_hist(idx,:)), '--', 'LineWidth', 1.0);

xlabel('Frame'); ylabel('|h|');

legend('True', 'LS'); title(sprintf('Coefficient %d magnitude', idx));

grid on;

end

end

-----

## 8\. MIMO Channel Capacity and Water-Filling Power Allocation

### Aim

[cite\_start]To simulate a MIMO channel and implement the *Water-Filling algorithm* to determine the *optimal power allocation* across the channel's orthogonal *eigenmodes, thereby **maximizing the achievable channel capacity*[cite: 2216].

### Key Observations

* [cite\_start]*Power Distribution:* The allocated power ($p_i$) is *highest for the strongest eigenmodes* and *zero for the weakest ones* that fall below the determined water level[cite: 2406, 2407].

* [cite\_start]*Capacity Maximization:* Capacity achieved with water-filling is the *maximum possible* due to the optimal power distribution based on *Channel State Information at the Transmitter (CSIT)*[cite: 2408, 2466].

### MATLAB Code (Listing 8.1 - Placeholder/Conceptual Logic)

matlab

% mimo_waterfill_example.m

% Example script to demonstrate MIMO Water-Filling power allocation

clear; clc; close all;

rng(2); % reproducible

%%% --System Parameters

Nr = 4; Nt = 4; % # receive and transmit antennas

Ptot = 10; % total transmit power

sigma2 = 1; % noise variance

% Generate a random Rayleigh MIMO channel matrix H ~ CN(0,1)

H = (randn(Nr, Nt) + 1j*randn(Nr, Nt)) / sqrt(2);

%% Simulation and Comparison

opts = struct();

opts.plot = true; % Option to plot the water-filling solution

opts.tol = 1e-9; % Tolerance for numerical methods

fprintf('Running k-search (closed-form) ---\n');

opts.method = 'ksearch';

% *** Placeholder function call for demonstration ***

out1 = run_waterfill_placeholder(H, Ptot, sigma2, opts);

print_results(out1);

fprintf('\n--- Running bisection (numeric) ---\n');

opts.method = 'bisection';

% *** Placeholder function call for demonstration ***

out2 = run_waterfill_placeholder(H, Ptot, sigma2, opts);

print_results(out2);

%% Helper Functions

% Small helper to print results nicely

function print_results(out)

fprintf('Method: %s\n', out.method);

fprintf('Status: %s\n', out.status);

fprintf('Eigenvalues (lambda):\n'); disp(out.lambdas.');

fprintf('Allocated powers p (per eigenmode): \n'); disp(out.p.');

fprintf('Water level mu = %.6f\n', out.mu);

fprintf('Sum allocated power = %.6f (Target Ptot= %.6f)\n', sum(out.p), out.Ptot);

fprintf('Capacity = %.6f bits/s/Hz\n', out.CbpsHz);

end

% Placeholder function that mimics the output structure for the example.

function out = run_waterfill_placeholder(H, Ptot, sigma2, opts)

% 1. SVD & Eigenvalues

[~, S, ~] = svd(H);

sigmas = diag(S);

lambdas = sigmas.^2;

lambdas_sorted = sort(lambdas, 'descend');

% 2. Simple Water-Filling calculation for demonstration (assuming all modes used).

% NOTE: This is a simplified calculation, full WF requires iterative k-search or bisection.

k_approx = length(lambdas_sorted); % assume all modes are used

mu_approx = (Ptot + sum(sigma2 ./ lambdas_sorted(1:k_approx))) / k_approx;

p_approx = max(0, mu_approx - sigma2 ./ lambdas_sorted);

% 3. Capacity

CbpsHz_approx = sum(log2(1 + lambdas_sorted .* p_approx / sigma2));

% 4. Output Structure

out.method = opts.method;

out.status = 'Success (Simulated)';

out.lambdas = lambdas_sorted;

out.p = p_approx;

out.mu = mu_approx;

out.CbpsHz = CbpsHz_approx;

out.Ptot = Ptot;

if opts.plot

% Generate a conceptual plot of the water-filling level

figure('Name', [opts.method ' Water-Filling'], 'NumberTitle', 'off');

% The concept is to plot allocated power (p_approx) as bars over the noise threshold

% The water level mu is related to the threshold by: Water Level = mu - sigma2/lambda_i

% The height of the bar is p_i, and the total conceptual 'height' is 1/lambda_i + p_i = mu

% The horizontal line plot is generally conceptual. Here we plot mu - sigma2/max(lambda) as a reference line.

bar(1:length(lambdas_sorted), p_approx, 'FaceColor', [0.7 0.7 1.0]); hold on;

% Conceptual noise line to show the 'water level' constraint

% The line should technically be at height mu, but often drawn at the active modes' level.

water_level_mu = mu_approx;

plot(xlim, [water_level_mu water_level_mu], '--r', 'LineWidth', 1.5);

title(['Water-Filling Power Allocation: ' opts.method]);

xlabel('Eigenmode Index (Strongest to Weakest)');

ylabel('Power Allocation ($p_i$)');

legend('Allocated Power $p_i$', 'Water Level $\mu$ (conceptual)');

grid on;

end

end

-----

## 9\. Orthogonal Frequency-Division Multiplexing (OFDM) Simulation

### Aim

[cite\_start]To simulate the fundamental working principle of an *Orthogonal Frequency-Division Multiplexing (OFDM) communication system, including **IFFT, Cyclic Prefix (CP) insertion, CP removal, **FFT*, and demodulation[cite: 2470].

### Key Observations

* [cite\_start]*Transmitter/Receiver Core:* The simulation confirms the effective use of the *IFFT/FFT pair* to multiplex data onto orthogonal subcarriers and the role of the *Cyclic Prefix (CP)* in the signal structure[cite: 2667].

* [cite\_start]*Received Constellation:* The symbols after the FFT cluster around the ideal *QPSK constellation* points, with scattering proportional to the AWGN[cite: 2639, 2662, 2663].

### MATLAB Code (Listing 9.1)

matlab

% ofdm_basic_sim.m

% Basic script to simulate a QPSK-OFDM system in an AWGN channel.

clear; clc; close all;

%%% --System Parameters

N = 64; % Number of subcarriers (FFT size)

cp_len = 16; % Cyclic Prefix length

M = 4; % Modulation order (QPSK)

numSymbols = 2000; % Number of OFDM symbols to simulate

bitsPerSym = log2(M);

numBits = N * numSymbols * bitsPerSym;

% Eb/NO in dB for AWGN channel

EbN0dB = 10; % Example value

%% 1. Transmitter

% Generate random bits

data = randi([0 1], numBits, 1);

% M-ary Modulation (QPSK)

txSym = pskmod(data, M, pi/4, 'gray', 'InputType', 'bit');

% Serial-to-Parallel conversion

txSymMatrix = reshape(txSym, N, numSymbols);

% IFFT (convert from frequency to time domain)

txIFFT = ifft(txSymMatrix, N);

% CP Insertion

txCP = [txIFFT(end - cp_len + 1 : end, :); txIFFT];

% Parallel-to-Serial conversion (OFDM signal)

txSignal = txCP(:);

%% 2. Channel

% Additive White Gaussian Noise (AWGN) channel

rxSignal = awgn(txSignal, EbN0dB, 'measured');

%% 3. Receiver

% Serial-to-Parallel conversion

rxSerial = reshape(rxSignal, N + cp_len, numSymbols);

% CP Removal

rxNoCP = rxSerial(cp_len + 1 : end, :);

% FFT (convert back to frequency domain)

rxFFT = fft(rxNoCP, N);

% Parallel-to-Serial conversion

rxSym = rxFFT(:);

% M-ary Demodulation (QPSK)

rxBits_dec = pskdemod(rxSym, M, pi/4, 'gray', 'OutputType', 'bit');

%% 4. Performance Metrics

% Calculate Bit Error Rate (BER)

[numErrors, BER] = biterr(data, rxBits_dec);

fprintf('\n- OFDM Simulation Results ---\n');

fprintf('Simulated Bits: %d\n', numBits);

fprintf('Total Bit Errors: %d\n', numErrors);

fprintf('Bit Error Rate (BER): %.6f\n', BER);

%% 5. Visualization (Constellation)

% Plot the received QPSK constellation after FFT

figure('Name', 'Received OFDM Constellation', 'NumberTitle', 'off');

scatter(real(rxSym(1:1000)), imag(rxSym(1:1000)), 'filled');

axis equal; grid on;

title(['Received QPSK Constellation (AWGN, Eb/NO = ' num2str(EbN0dB) ' dB)']);

xlabel('In-Phase'); ylabel('Quadrature');

-----

## 10\. Orthogonal Space-Time Block Codes (OSTBC) and Alamouti Code

### Aim

[cite\_start]To implement and simulate the performance of *Orthogonal Space-Time Block Codes (OSTBC), specifically the **Alamouti Code ($2\times1$ configuration), to demonstrate the principles of **transmit diversity* and analyze *BER* performance in a Rayleigh fading channel[cite: 2673].

### Key Observations

* [cite\_start]*Diversity Gain:* The Alamouti code demonstrates a *2nd order diversity gain* (for $N_t=2, N_r=1$), showing a *much steeper BER curve* than a single-antenna system[cite: 2964].

* [cite\_start]*Code Functionality:* The manual implementation confirms the principle of *orthogonal encoding* and the use of *Maximum Ratio Combining (MRC)-like linear decoding* for optimal detection[cite: 2698, 2964].

### MATLAB Code (Listing 10.1 - Alamouti 2x1 Manual Implementation)

matlab

% alamouti_sim.m Alamouti 2x1 STBC BER simulation

clear; close all; rng(1);

M = 4; k = log2(M); numBits = 2e5;

bits = randi([0 1], numBits, 1);

% QPSK modulation (symbols normalized to unit average energy)

sym = (2*bits(1:2:end)-1) + 1j*(2*bits(2:2:end)-1);

sym = sym / sqrt(2);

EbN0dB = 0:2:18;

ber = zeros(size(EbN0dB));

for idx = 1:length(EbN0dB)

EbNO = 10^(EbN0dB(idx)/10);

EsNO = EbNO * k;

noiseVar = 1 / (2 * EsNO); % Noise variance per real/imaginary dimension

rxBits = zeros(numBits, 1);

recvPtr = 1;

for n = 1:2:length(sym)-1

s1 = sym(n);

s2 = sym(n+1);

% Generate Rayleigh channel coefficients (h1, h2)

h = (randn(1, 2) + 1j*randn(1, 2)) / sqrt(2);

% Generate complex AWGN noise

n1 = sqrt(noiseVar) * (randn + 1j*randn);

n2 = sqrt(noiseVar) * (randn + 1j*randn);

% Received Signals (Two time slots)

% Time 1: y1 = h1*s1 + h2*s2 + n1

y1 = h(1) * s1 + h(2) * s2 + n1;

% Time 2: y2 = h1*(-s2*) + h2*(s1*) + n2

y2 = h(1) * (-conj(s2)) + h(2) * conj(s1) + n2;

% Maximal Ratio Combining (MRC)-like Decoding

denom = abs(h(1))^2 + abs(h(2))^2;

% s1_hat = (h1*y1 + h2*y2*) / G

s1_hat = (conj(h(1)) * y1 + h(2) * conj(y2)) / denom;

% s2_hat = (h2*y1 - h1*y2*) / G

s2_hat = (conj(h(2)) * y1 - h(1) * conj(y2)) / denom;

% Hard Decision (Real and Imaginary parts)

rxBits(recvPtr) = real(s1_hat) > 0;

rxBits(recvPtr+1) = imag(s1_hat) > 0;

rxBits(recvPtr+2) = real(s2_hat) > 0;

rxBits(recvPtr+3) = imag(s2_hat) > 0;

recvPtr = recvPtr + 4;

end

ber(idx) = mean(rxBits(1:numBits) ~= bits);

fprintf('Eb/NO %d dB, BER = %g\n', EbN0dB(idx), ber(idx));

end

figure;

semilogy(EbN0dB, ber, 'o-', 'LineWidth', 1.4);

xlabel('E_b/N_0 (dB)'); ylabel('BER'); grid on;

title('Alamouti 2x1 STBC QPSK over Rayleigh Fading');