Telecommunication MATLAB Projects (20+ Ideas with Complete Code)

% 1. Synthesize Telemetry Sensor Data & Spreading Code
telemetry_data = randi([0 1], 100, 1); % 100-bit telemetry frame
pnGen = comm.PNSequence('Polynomial', [1 0 0 0 0 1 1], 'SamplesPerFrame', 3100);
chip_seq = pnGen(); % 31 chips per bit processing gain

% 2. Direct Sequence Spread Spectrum (DSSS) Modulation
bpskMod = comm.BPSKModulator();
tx_symbols = bpskMod(telemetry_data);
tx_spread = repelem(tx_symbols, 31) .* (1 - 2*chip_seq);

% 3. Transmit through Multipath AWGN Channel
rx_spread = awgn(tx_spread, 5, 'measured');

% 4. Correlation Despreading & Demodulation
despread_sig = rx_spread .* (1 - 2*chip_seq);
rx_integrated = sum(reshape(despread_sig, 31, []), 1)';
rx_bits = rx_integrated > 0;
[numErrors, ber] = biterr(telemetry_data, rx_bits);
fprintf('Telemetry Link BER: %.4e | Total Bit Errors: %d\n', ber, numErrors);

N_sc = 64; N_cp = 16; M = 16; % 16-QAM, 64 subcarriers
data_bits = randi([0 M-1], N_sc, 100);
mod_syms = qammod(data_bits, M, 'UnitAveragePower', true);

% IFFT Modulation & Cyclic Prefix Insertion
ifft_out = ifft(mod_syms, N_sc, 1);
tx_ofdm = [ifft_out(end-N_cp+1:end, :); ifft_out];

% Multipath Channel & AWGN
h = [0.8; 0.4; 0.2; 0.1]; % 4-tap channel impulse response
rx_ofdm = filter(h, 1, tx_ofdm(:));
rx_noisy = awgn(rx_ofdm, 20, 'measured');
rx_matrix = reshape(rx_noisy(1:end-mod(length(rx_noisy), N_sc+N_cp)), N_sc+N_cp, []);

% Remove CP & FFT Demodulation
rx_no_cp = rx_matrix(N_cp+1:end, :);
fft_out = fft(rx_no_cp, N_sc, 1);
H_freq = fft(h, N_sc);
equalized_syms = fft_out ./ H_freq; % Zero-Forcing Equalizer

scatterplot(equalized_syms(:)); title('Equalized 16-QAM OFDM Constellation');

% 1. Parameters
N_bits = 1000; L_code = 63; % 63-chip m-sequence code
data = randi([0 1], N_bits, 1);
bpsk_data = 2*data - 1;

% 2. Generate 63-chip PN Spreading Sequence
pn = comm.PNSequence('Polynomial', 'x^6 + x^1 + 1', 'SamplesPerFrame', N_bits*L_code);
pn_seq = 2*pn() - 1;
spread_sig = repelem(bpsk_data, L_code) .* pn_seq;

% 3. Add Strong Narrowband Jammer Tone
t = (1:length(spread_sig))';
jammer = 2.5 * cos(2*pi*0.05*t); % High-power tone
rx_sig = awgn(spread_sig + jammer, 10, 'measured');

% 4. Despreading & Matched Filtering
despread = rx_sig .* pn_seq;
recovered = sum(reshape(despread, L_code, N_bits), 1)' > 0;
[~, ber] = biterr(data, recovered);
fprintf('DSSS Anti-Jamming BER: %.4f (Processing Gain: %.1f dB)\n', ber, 10*log10(L_code));

% Pulse Shaping & Eye Diagram Educational Demo
M = 4; sps = 8; rolloff = 0.35;
tx_filter = comm.RaisedCosineTransmitFilter('RolloffFactor', rolloff, ...
    'FilterSpanInSymbols', 6, 'OutputSamplesPerSymbol', sps);

symbols = qammod(randi([0 M-1], 1000, 1), M, 'UnitAveragePower', true);
tx_pulses = tx_filter(symbols);
rx_noisy = awgn(tx_pulses, 18, 'measured');

% Render Interactive Eye Diagram
eyediagram(rx_noisy(1:1000), 2*sps, 2, sps/2);
title(['Eye Diagram (QPSK, Rolloff = ', num2str(rolloff), ', SNR = 18 dB)']);

% 1. SAR Parameters & Fast/Slow Time Grid
c = 3e8; fc = 5e9; B = 100e6; fs = 200e6;
N_range = 512; N_azimuth = 256;
K_r = B / (N_range/fs); % Chirp rate

% 2. Synthesize Raw SAR Echo from Point Targets
raw_echo = randn(N_range, N_azimuth) * 0.05 + 1j*randn(N_range, N_azimuth)*0.05;
raw_echo(250:260, 120:130) = raw_echo(250:260, 120:130) + 10; % Target 1

% 3. Range Compression via 1D FFT Matched Filter
ref_chirp = exp(-1j * pi * K_r * ((0:N_range-1)/fs).^2).';
range_compressed = ifft(fft(raw_echo, [], 1) .* conj(fft(ref_chirp, N_range, 1)), [], 1);

% 4. Azimuth FFT & Matched Filtering
sar_image = abs(fftshift(fft(range_compressed, [], 2), 2));
figure; imagesc(20*log10(sar_image)); colormap('jet'); colorbar;
title('Reconstructed 2D SAR Radar Image (dB)'); xlabel('Azimuth Bin'); ylabel('Range Bin');

% 1. Generate Smart Meter Power Load Profile (kW)
hours = 1:24;
load_profile = 1.5 + 2.0*sin((hours-6)*pi/12) + 0.3*randn(size(hours));
load_profile(load_profile < 0.2) = 0.2;

% 2. Differential Privacy Laplace Noise Injection
epsilon = 0.8; % Privacy budget
sensitivity = max(diff(load_profile));
b = sensitivity / epsilon; % Laplace scale parameter
u = rand(size(hours)) - 0.5;
laplace_noise = -b * sign(u) .* log(1 - 2*abs(u));

obfuscated_load = load_profile + laplace_noise;

figure; plot(hours, load_profile, 'b-o', hours, obfuscated_load, 'r--s', 'LineWidth', 1.5);
legend('True Private Meter Load', 'Obfuscated Telemetry (\epsilon=0.8)');
xlabel('Hour of Day'); ylabel('Power Demand (kW)'); grid on;

% 1. Grid & Soliton Parameters
N = 2048; dt = 10e-15; T = (-N/2:N/2-1)*dt;
T0 = 50e-15; P0 = 1000; % 50 fs pulse, 1 kW peak
gamma = 0.05; beta2 = -20e-27; % Anomalous dispersion

% 2. Initial Sech Pulse Field
A = sqrt(P0) * sech(T / T0);
w = 2*pi*(-N/2:N/2-1)/(N*dt);
dz = 0.001; z_max = 0.1; % 10 cm fiber

% 3. Split-Step Fourier Method (Linear Dispersion + Nonlinear SPM)
D_op = exp(1j * 0.5 * beta2 * (w.^2) * dz);
for z = 0:dz:z_max
    A = A .* exp(1j * gamma * abs(A).^2 * (dz/2));
    A = ifft(fftshift(D_op .* fftshift(fft(A))));
    A = A .* exp(1j * gamma * abs(A).^2 * (dz/2));
end

figure; plot(w/(2*pi*1e12), 10*log10(abs(fftshift(fft(A))).^2));
xlabel('Frequency Offset (THz)'); ylabel('Spectral Power (dB)'); title('Supercontinuum Broadening');

% 1. Synthesize 3D Urban Space Coordinates
[X, Y, Z] = meshgrid(-50:5:50, -50:5:50, 0:5:30);
pts = [X(:), Y(:), Z(:)];

% 2. Compute 3D Antenna Radiation Field (Directivity & Path Loss)
dist = sqrt(pts(:,1).^2 + pts(:,2).^2 + pts(:,3).^2) + 1;
theta = atan2(sqrt(pts(:,1).^2 + pts(:,2).^2), pts(:,3));
radiation_pattern = (cos(theta).^2) ./ (dist.^1.5);
rssi_values = 10*log10(radiation_pattern / max(radiation_pattern)) - 40;

% 3. Render 3D RF Field Intensity PointCloud
ptCloud = pointCloud(pts, 'Intensity', single(rssi_values));
figure; pcshow(ptCloud); colormap('turbo'); colorbar;
title('3D Spatial RF Signal Intensity PointCloud (dBm)');

% 1. Define BPM-to-Corrector Magnet Response Matrix R
N_bpm = 32; N_corr = 16;
R = randn(N_bpm, N_corr); % Measured orbit response matrix

% 2. Singular Value Decomposition (SVD) Orbit Inversion
[U, S, V] = svd(R, 'econ');
s_diag = diag(S);
k_modes = 10; % Retain top 10 dominant spatial modes
R_inv = V(:,1:k_modes) * diag(1./s_diag(1:k_modes)) * U(:,1:k_modes)';

% 3. Simulate Orbit Perturbation & Closed-Loop Feedback
measured_orbit = 0.5 * sin((1:N_bpm)' * 2*pi/8) + 0.05*randn(N_bpm, 1);
corrector_currents = -R_inv * measured_orbit;
corrected_orbit = measured_orbit + R * corrector_currents;

fprintf('Orbit Distortion RMS Before: %.3f mm | After SVD Correction: %.3f mm\n', ...
    rms(measured_orbit), rms(corrected_orbit));

% 1. Synthesize 2nd-Derivative Gaussian Monocycle Pulse
fc = 4e9; fs = 50e9; tc = 0.5e-9;
t_pulse = -1e-9:1/fs:1e-9;
p = (1 - 4*pi*(t_pulse/tc).^2) .* exp(-2*pi*(t_pulse/tc).^2); % UWB Monocycle

% 2. Time-Hopping Pulse Position Modulation (TH-PPM)
bits = [1 0 1 1 0]; delta = 0.2e-9; % PPM delay
tx_signal = [];
for b = bits
    t_shift = b * delta;
    tx_signal = [tx_signal, zeros(1, 100), p]; % Pulse frame
end

% 3. Channel & Template Correlation Detection
rx_noisy = awgn(tx_signal, 10, 'measured');
corr_out = xcorr(rx_noisy, p);

figure; subplot(2,1,1); plot(tx_signal); title('Transmitted TH-PPM UWB Pulses');
subplot(2,1,2); plot(corr_out); title('Matched Filter Correlation Demodulation');

% 1. Configure 5G NR Carrier & PDSCH Configuration
carrier = nrCarrierConfig('NSizeGrid', 52, 'SubcarrierSpacing', 30);
pdsch = nrPDSCHConfig('Modulation', '64QAM', 'NumLayers', 2);
[pdschIndices, pdschInfo] = nrPDSCHIndices(carrier, pdsch);

% 2. Generate Information Bits & PDSCH Waveform
dataBits = randi([0 1], pdschInfo.G, 1);
pdschSymbols = nrPDSCH(carrier, pdsch, dataBits);

% 3. Map to 5G NR Resource Grid
grid = nrResourceGrid(carrier, 2);
grid(pdschIndices) = pdschSymbols;

% 4. CDL Multipath Channel Transmission
cdl = nrCDLChannel('DelayProfile', 'CDL-C', 'DelaySpread', 100e-9);
cdl.SampleRate = 30.72e6;
rxGrid = awgn(grid, 22, 'measured');

fprintf('5G NR PDSCH Grid Configured: 52 PRBs @ 30 kHz SCS, 64-QAM, Rate: %.2f Mbps\n', ...
    length(dataBits)/(1e-3 * 1e6));

% 1. Configure 64-Antenna Base Station Array & 4 Single-Antenna Users
Nt = 64; K = 4;
H = (randn(K, Nt) + 1j*randn(K, Nt)) / sqrt(2); % Rayleigh Channel Matrix

% 2. Multi-User Zero-Forcing (ZF) Digital Precoder
W_zf = H' / (H * H');
W_zf = W_zf ./ sqrt(sum(abs(W_zf).^2, 1)); % Power normalization

% 3. Compute Multi-User Sum-Rate Capacity
SNR_dB = 0:5:30;
sum_rate = zeros(size(SNR_dB));
for i = 1:length(SNR_dB)
    snr = 10^(SNR_dB(i)/10);
    signal_power = abs(diag(H * W_zf)).^2 * snr;
    sinr = signal_power ./ (1 + sum(abs(H*W_zf).^2, 2) - signal_power);
    sum_rate(i) = sum(log2(1 + sinr));
end

figure; plot(SNR_dB, sum_rate, 'b-o', 'LineWidth', 2);
xlabel('SNR (dB)'); ylabel('Sum Spectral Efficiency (bps/Hz)');
title('64x4 Massive MIMO Zero-Forcing Sum Capacity'); grid on;

% 1. UE Movement Trajectory between BS1 (0m) and BS2 (1000m)
d = 10:5:990; % User position along baseline
PL1 = 128.1 + 37.6*log10(d/1000) + randn(size(d))*3.5;
PL2 = 128.1 + 37.6*log10((1000-d)/1000) + randn(size(d))*3.5;
RSRP1 = 43 - PL1; RSRP2 = 43 - PL2;

% 2. Event A3 Evaluation with Hysteresis & Time-to-Trigger
Hyst = 3.0; TTT_steps = 5;
ho_triggered = false; serving_cell = ones(size(d));
counter = 0;

for k = 1:length(d)
    if (RSRP2(k) - RSRP1(k)) > Hyst
        counter = counter + 1;
        if counter >= TTT_steps && ~ho_triggered
            serving_cell(k:end) = 2;
            ho_triggered = true;
        end
    else
        counter = 0;
    end
end

figure; plot(d, RSRP1, 'b', d, RSRP2, 'r', 'LineWidth', 1.5);
hold on; stairs(d, serving_cell*10 - 100, 'k--', 'LineWidth', 2);
legend('BS1 RSRP', 'BS2 RSRP', 'Serving Cell ID'); xlabel('Distance (m)'); ylabel('RSRP (dBm)');

% 1. Ka-Band LEO Satellite Link Budget Calculation
fc = 20e9; Pt_dBW = 15; Gt_dBi = 32; Gr_dBi = 42; % Tx/Rx parameters
h_orbit = 600e3; c = 3e8;
lambda = c / fc;
FSPL_dB = 20*log10(4*pi*h_orbit/lambda); % Free space path loss
Atm_loss_dB = 1.5; Rain_loss_dB = 3.0;

EIRP = Pt_dBW + Gt_dBi;
Pr_dBW = EIRP + Gr_dBi - FSPL_dB - Atm_loss_dB - Rain_loss_dB;
fprintf('Satellite Link Margin: EIRP = %.1f dBW | Received Power = %.2f dBm\n', ...
    EIRP, Pr_dBW + 30);

% 2. Dynamic LEO Doppler Shift Synthesis
t = 0:0.001:5; v_rel = 7500; % 7.5 km/s orbital speed
f_doppler = (v_rel / c) * fc * (1 - (2*t/5)); % Linear Doppler profile
carrier_sig = exp(1j * 2*pi * cumsum(f_doppler)*0.001);

% 1. Optical WDM Transmitter (8 Channels @ 10 Gbps)
num_ch = 8; bit_rate = 10e9; N_bits = 2048;
data = randi([0 1], N_bits, num_ch);
opt_pulse = repelem(data, 16); % NRZ electrical drive signal

% 2. Fiber Channel: 80 km Span with Attenuation & CD
alpha_dB = 0.2; L_span = 80; % km
D_cd = 16; % ps/(nm*km) dispersion parameter
attenuation = 10^(-(alpha_dB * L_span)/20);

% 3. EDFA Optical Amplifier Simulation
edfa_gain_dB = alpha_dB * L_span;
edfa_gain = 10^(edfa_gain_dB/20);
noise_fig_dB = 5;
ase_noise = randn(size(opt_pulse)) * 0.02;

rx_optical = opt_pulse * attenuation * edfa_gain + ase_noise;
eyediagram(rx_optical(1:500, 1), 32); title('Channel 1 Optical Eye Diagram (80 km Post-EDFA)');

% 1. Simulation Parameters
N_samples = 1000; num_trials = 5000; snr_dB = -8;
snr = 10^(snr_dB/10); Pfa_target = 0.01:0.02:0.99;

% 2. Theoretical Neyman-Pearson Decision Threshold
sigma_w2 = 1;
gamma_th = sigma_w2 * (qfuncinv(Pfa_target)/sqrt(N_samples) + 1);

% 3. Monte Carlo Sensing Simulation
Pd_sim = zeros(size(Pfa_target));
for k = 1:length(Pfa_target)
    noise_only = sqrt(sigma_w2/2) * (randn(N_samples, num_trials) + 1j*randn(N_samples, num_trials));
    pu_signal = sqrt(snr/2) * (randn(N_samples, num_trials) + 1j*randn(N_samples, num_trials));
    test_statistic = sum(abs(pu_signal + noise_only).^2, 1) / N_samples;
    Pd_sim(k) = sum(test_statistic > gamma_th(k)) / num_trials;
end

figure; plot(Pfa_target, Pd_sim, 'b-o', 'LineWidth', 2); grid on;
xlabel('Probability of False Alarm (P_{fa})'); ylabel('Probability of Detection (P_d)');
title(['Cognitive Radio ROC Curve (SNR = ', num2str(snr_dB), ' dB)']);

% 1. 5G NR Polar Code Parameters (K=64 message bits, E=128 encoded bits)
K = 64; E = 128; L = 8; % SCL List Length
msg_bits = randi([0 1], K, 1);

% 2. 5G Polar Encoding
encoded_bits = nrPolarEncode(msg_bits, E);
bpsk_symbols = 1 - 2*encoded_bits;

% 3. AWGN Channel & Log-Likelihood Ratio (LLR)
EbNo = 4; % dB
snr = EbNo + 10*log10(K/E) + 10*log10(2);
rx_noisy = awgn(bpsk_symbols, snr, 'measured');
llr = 2 * rx_noisy * 10^(snr/10);

% 4. Successive Cancellation List (SCL) Decoding
decoded_bits = nrPolarDecode(llr, K, E, L);
[numErr, ber] = biterr(msg_bits, decoded_bits);
fprintf('5G Polar Code Decoding: Errors = %d | BER = %.2e (SCL List = %d)\n', numErr, ber, L);

% 1. Gamma-Gamma Atmospheric Turbulence Parameters
alpha = 4.2; beta = 1.4; % Moderate turbulence parameters
I_samples = gamrnd(alpha, 1/alpha, [10000, 1]) .* gamrnd(beta, 1/beta, [10000, 1]);

% 2. OOK Modulated Signal Transmission
data_bits = randi([0 1], 10000, 1);
tx_laser = data_bits .* I_samples; % Scintillation modulated intensity

% 3. Photodetector & AWGN Thermal Noise
rx_signal = awgn(tx_laser, 15, 'measured');
decision_threshold = 0.5 * mean(I_samples);
detected_bits = rx_signal > decision_threshold;

[~, ber] = biterr(data_bits, detected_bits);
fprintf('FSO Link BER under Gamma-Gamma Scintillation: %.4e\n', ber);

% 1. LoRa Parameters (SF = 7, Bandwidth = 125 kHz)
SF = 7; BW = 125e3; Fs = BW; N = 2^SF;
symbol = 42; % Data symbol (0 to N-1)

% 2. Generate Cyclically Shifted Chirp Waveform
k = 0:N-1;
f_inst = mod(symbol + k, N) * (BW/N);
phi = 2*pi * cumsum(f_inst)/Fs;
tx_chirp = exp(1j * phi);

% 3. AWGN Channel (-12 dB Deep Sub-Noise SNR)
rx_noisy = awgn(tx_chirp, -12, 'measured');

% 4. De-Chirping Conjugate Multiplication & FFT Demodulation
base_downchirp = exp(-1j * 2*pi * (0.5*BW/N * (k.^2)) / Fs);
dechirped = rx_noisy .* base_downchirp;
spec = abs(fft(dechirped, N));
[~, detected_symbol] = max(spec);
detected_symbol = detected_symbol - 1;

fprintf('LoRa CSS Demodulation: Transmitted = %d | Detected = %d (SNR: -12 dB)\n', ...
    symbol, detected_symbol);

% 1. Generate QPSK Symbols with Barker Preamble
preamble = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1]'; % 13-bit Barker
payload = randi([0 3], 200, 1);
payload_qpsk = pskmod(payload, 4, pi/4);
tx_frame = [preamble; payload_qpsk];

% 2. Channel Impairments: Carrier Frequency Offset (CFO) & Noise
cfo = 0.02; % Normalized frequency offset
t = (1:length(tx_frame))';
rx_impaired = awgn(tx_frame .* exp(1j*2*pi*cfo*t), 18, 'measured');

% 3. Closed-Loop Carrier Phase & Frequency Synchronizer (Costas Loop)
carrierSync = comm.CarrierSynchronizer('SamplesPerSymbol', 1, ...
    'Modulation', 'QPSK', 'DampingFactor', 0.707, 'NormalizedLoopBandwidth', 0.01);
[rx_synced, phase_err] = carrierSync(rx_impaired);

scatterplot(rx_synced(50:end)); title('SDR Synchronized QPSK Constellation');