hfc_setup_1.m

% real-life test using PlutoSDR

clear all;
clc;
close all;

M = 32;  % QAM orders (up to 1024) 
k = log2(M);
N_fft = 512;  %OFDM block dimension
N_cp = 16;  %Cyclic prefix length (varying from 8 to 32)
Fs = 12.8e6; %sample rate
% %Fs = 50e6; %sample rate
% %Fs = 40e6; %sample rate
% %Fs = 2e6; %sample rate USB2

guard_l = 18; %number if guard subcarriers on left side
guard_r = 17; %on right side

%% TX construction
N_data=N_fft-guard_l-guard_r-1; %number of active data subcarriers
nb_blocks=1000;

N_symbol=N_data*nb_blocks;

tx_bits = randi([0 1], N_symbol * k, 1);
tx_symb = qammod(tx_bits, M, InputType='bit', UnitAveragePower=true);

%organization in blocks
data_matrix=reshape(tx_symb, N_data, nb_blocks); %N_data lines et nb_blocks columns
spectrum_centered=zeros(N_fft, nb_blocks); %creation spectrum
%separation of symbols for negative and positive frequencies
half_active=N_data/2;
spectrum_centered(guard_l+1 : N_fft/2, :) = data_matrix(1:half_active, :);
spectrum_centered(N_fft/2+2 : N_fft-guard_r, :) = data_matrix(half_active+1:end, :);

% Passage temporel + CP pour les données
tx_signal = ifft(ifftshift(spectrum_centered, 1), N_fft);
cp_samples = tx_signal(end-N_cp+1:end, :);
tx_ofdm_cp = [cp_samples; tx_signal];

% Preamble reneration
preamble_bits=randi([0 1],N_data* k, 1);
preamble_symb=qammod(preamble_bits,M, InputType='bit',UnitAveragePower=true);
spectrum_preamble = zeros(N_fft, 1);
spectrum_preamble(guard_l+1 : N_fft/2) = preamble_symb(1 : half_active);
spectrum_preamble(N_fft/2+2 : N_fft-guard_r) = preamble_symb(half_active+1 : end);

tx_preamble_time = ifft(ifftshift(spectrum_preamble, 1), N_fft);
cp_preamble = tx_preamble_time(end-N_cp+1 : end);
tx_preamble_block = [cp_preamble; tx_preamble_time]; % block with cp 

% construction of the complete package
tx_packet_matrix = [tx_preamble_block, tx_ofdm_cp];
tx_packet_serial = tx_packet_matrix(:);

%% PlutoSDR setup
% need a `buffer_tx` vector in workspace
buffer_tx = tx_packet_serial; % which ? should be tx_packet_serial

% TODO: value ?
Rsamp = Fs; % base sample rate
%tx_gain_dB = -10;
%rx_gain_dB = 10;
%center_freq = 1e9;
%while 1==1
% --- Transmitter system object ---
txPluto = sdrtx('Pluto',...
    'RadioID', 'usb:0',...
    'Gain', -10,... % -90 to 0 dB
    'CenterFrequency', 820e6,... % 335e6 to 3.8e9 [Hz]
    'BasebandSampleRate', Rsamp); % 60e3 to 60e6 [Hz]

% --- Receiver system object ---
rxPluto = sdrrx('Pluto',...
    'RadioID', 'usb:0',...
    'CenterFrequency', 820e6,... % 335e6 to 3.8e9 [Hz]
    'GainSource', 'Manual',... % Manual, AGC Slow Attack
    'Gain', 10,... % -4 to 71 [dB]
    'BasebandSampleRate', Rsamp,... % 60e3 to 60e6 [Hz]
    'EnableBurstMode', true,...
    'NumFramesInBurst', 1,...
    'SamplesPerFrame', 3*size(tx_packet_serial,1),...%200000,... % <-------------
    'OutputDataType', 'double' );

%% PlutoSDR transmit

% --- Periodic transmission ---
txPluto.transmitRepeat(buffer_tx);

% ! Note that since the buffer is periodically transmitted, the transients 
% effects at the beginning and at the end of the signal due to the filter 
% memory should first be removed. !
figure;
%% PlutoSDR receive

% --- Burst reception ---
% Since the transmitted signal is infinitely periodic, no difficult sample 
% synchronisation needs to be foreseen when recording the signal.

% https://nl.mathworks.com/help/comm/plutoradio/ref/comm.sdrrxpluto-system-object.html#d127e7772
[data_rx,datavalid,overflow] = rxPluto();

if datavalid
    disp('Received sample is valid');
else
    disp('Received sample is not valid');
end

if (overflow)
    disp('Samples dropped');
end

% data_rx: Received data, returned as a vector of complex ('double') values

%% RX processing

% Discard the first and last few thousand samples to bypass transients
N_discard = 5000; 
data_rx = data_rx(N_discard+1 : end-N_discard);

rx_packet_serial = data_rx; % this is our 'rx_signal_noisy' from our simulation code
% TODO: process rx_packet_serial
% ... cross-corr ... sync ... channel estimation ... raw zero-forcing ...
% denoising ... symbol reconstruction ... qamdemod ... done (compare with tx)

rx_signal_noisy = rx_packet_serial;
% cross-correlation
[C, lags] = xcorr(rx_signal_noisy, tx_preamble_time); % FIXME: should we xcorr with the preamble block (including its CP) ?!
%[C, lags] = xcorr(rx_signal_noisy, tx_preamble_block);
valid_idx = lags >= 0;
C_fixed = C(valid_idx);
lags_fixed = lags(valid_idx);

[~, max_idx] = max(abs(C_fixed));
estimated_ToA = lags_fixed(max_idx)-N_cp -8; % TODO: less/more than N_cp ?

% prevent negative indexing if the preamble arrived too early in the buffer
if estimated_ToA < 0
    estimated_ToA = 0; % So that (estimated_ToA + 1) safely targets index 1
    disp('Warning: Preamble arrived very early, ToA clamped to start of buffer.');
end

% synchronization
rx_synchronized = rx_signal_noisy(estimated_ToA + 1 : end);

L_bloc=N_fft+N_cp;
N_total_samples =L_bloc * (nb_blocks + 1);

% prevent reshape errors if the capture ran out of samples at the end
if length(rx_synchronized) < N_total_samples
    error('Not enough samples remaining to reshape! Increase SamplesPerFrame in rxPluto. Or re-run script.');
end

rx_packet_matrix_est = reshape(rx_synchronized(1:N_total_samples), L_bloc, nb_blocks + 1);

%% channel estimation
% extraction of the received preamble
rx_preamble_signal = rx_packet_matrix_est(N_cp+1:end, 1);
rx_preamble_fft = fftshift(fft(rx_preamble_signal, N_fft), 1);

active_idx = [guard_l+1 : N_fft/2, N_fft/2+2 : N_fft-guard_r];

% raw zero-forcing
H_raw_active = rx_preamble_fft(active_idx) ./ spectrum_preamble(active_idx); % $\underline{\hat{H}}_{raw} = \underline{\underline{\Lambda}}_I^{-1} \cdot \underline{x}$

% % ---------
% % Reconstruct a full 512-point raw frequency channel (pad inactive subcarriers with 0)
% H_raw_full = zeros(N_fft, 1);
% H_raw_full(active_idx) = H_raw_active;
% 
% % Convert to time domain to see the raw impulse response
% h_time_full = ifft(H_raw_full);
% 
% figure('Name', 'Time Domain Channel Impulse Response');
% plot(abs(h_time_full), 'LineWidth', 1.5);
% title('Raw Channel Impulse Response');
% xlabel('Tap Index');
% ylabel('Magnitude');
% xlim([1 N_fft]);
% grid on;
% % ---------

% denoising
L_taps = N_cp*2;  % TODO: vary ?
% building the truncated Fourier matrix (slide 62):
[Columns, Rows] = meshgrid(0:L_taps-1, active_idx-1); % Rows are active freqs, Columns are L valid channel taps
F_active = exp(-1i * 2 * pi * Rows .* Columns / N_fft); % $\overline{\underline{\underline{F}}}_Q$ ($N_{active} \times L$ matrix)

h_estim_time = F_active \ H_raw_active; % $F_{active} \cdot h = H_{raw}$
% this does $\underline{\hat{h}} = (\overline{\underline{\underline{F}}}_Q^H \cdot \overline{\underline{\underline{F}}}_Q)^{-1} \cdot \overline{\underline{\underline{F}}}_Q^H \cdot \underline{\hat{H}}_{raw}$

H_denoised_active = F_active * h_estim_time; % back to FD ($\underline{\hat{H}} = \sqrt{Q} \underline{\underline{F}}_Q \underline{\hat{h}}$)

%% symbol reconstruction and qamdemod
rx_data_blocs = rx_packet_matrix_est(N_cp+1:end, 2:end);
rx_data_freq = fftshift(fft(rx_data_blocs, N_fft), 1);

rx_active = [rx_data_freq(guard_l+1 : N_fft/2, :); rx_data_freq(N_fft/2+2 : N_fft-guard_r, :)];

symb_no_eq = rx_active(:);
symb_raw = rx_active(:) ./ repmat(H_raw_active, nb_blocks, 1);
symb_denoised = rx_active(:) ./ repmat(H_denoised_active, nb_blocks, 1);

rx_bits_no_eq = qamdemod(symb_no_eq, M, OutputType='bit', UnitAveragePower=true);
rx_bits_raw_eq = qamdemod(symb_raw, M, OutputType='bit', UnitAveragePower=true);
rx_bits_denoised = qamdemod(symb_denoised, M, OutputType='bit', UnitAveragePower=true);

%%

% % Constellation 
% figure('Name', 'Constellation');
% scatter(real(symb_no_eq), imag(symb_no_eq), 10, 'b', 'filled'); hold on;
% title(sprintf('QAM-%d Constellation at Receiver', M));
% xlabel('In-Phase'); ylabel('Quadrature');
% legend('Received Symbols');
% grid on; axis square;
% % Constellation 
% figure('Name', 'Constellation');
% scatter(real(symb_raw), imag(symb_raw), 10, 'b', 'filled'); hold on;
% title(sprintf('QAM-%d Constellation at Receiver', M));
% xlabel('In-Phase'); ylabel('Quadrature');
% legend('Received Symbols');
% grid on; axis square;
% Constellation 
%figure('Name', 'Constellation');
scatter(real(symb_denoised), imag(symb_denoised), 10, 'b', 'filled'); hold on;
title(sprintf('QAM-%d Constellation at Receiver', M));
xlabel('In-Phase'); ylabel('Quadrature');
legend('Received Symbols');
grid on; axis square;
%end

%%

% compute SNR for each equalization method 
% TODO: maybe use symb instead? tx_symb and symb_no_eq etc
%SNR_no_eq = 10 * log10(var(tx_bits) / var(tx_bits - rx_bits_no_eq))
%SNR_raw_eq = 10 * log10(var(tx_bits) / var(tx_bits - rx_bits_raw_eq))
%SNR_denoised_eq = 10 * log10(var(tx_bits) / var(tx_bits - rx_bits_denoised))
SNR_eff_no_eq = 10 * log10(var(tx_symb) / var(tx_symb - symb_no_eq))
SNR_eff_raw_eq = 10 * log10(var(tx_symb) / var(tx_symb - symb_raw))
SNR_eff_denoised_eq = 10 * log10(var(tx_symb) / var(tx_symb - symb_denoised))

% bar chart for SNRs
%figure('Name','SNRs EQ methods real-life');
figure('Name','Real-life experiment data');
sgtitle('OFDM Link Performance for Different Equalization Methods'); % <--- Adds the main figure title
subplot(1,2,1);
b_snr = bar([SNR_eff_no_eq, SNR_eff_raw_eq, SNR_eff_denoised_eq]);
% 1. Enable individual bar coloring
b_snr.FaceColor = 'flat';
% 2. Assign RGB colors (using MATLAB's default nice color palette)
b_snr.CData(1,:) = [0.8500 0.3250 0.0980]; % Red/Orange for No EQ
b_snr.CData(2,:) = [0.0000 0.4470 0.7410]; % Blue for Raw EQ
b_snr.CData(3,:) = [0.4660 0.6740 0.1880]; % Green for Denoised EQ
%xlabel('Equalization Method');
ylabel('Effective SNR (dB)');
title('Effective SNR Comparison');
xticks(1:3);
xticklabels({'No EQ', 'Raw EQ', 'Denoised EQ'});
grid on;

% compute BER
BER_no_eq = sum(tx_bits ~= rx_bits_no_eq) / length(tx_bits)
BER_raw_eq = sum(tx_bits ~= rx_bits_raw_eq) / length(tx_bits)
BER_denoised_eq = sum(tx_bits ~= rx_bits_denoised) / length(tx_bits)
%BER_no_eq = sum(tx_symb ~= symb_no_eq) / length(tx_symb)
%BER_raw_eq = sum(tx_symb ~= symb_raw) / length(tx_symb)
%BER_denoised_eq = sum(tx_symb ~= symb_denoised) / length(tx_symb)

% bar chart for BERs
subplot(1,2,2);
%figure('Name','BERs EQ methods real-life');
b_ber = bar([BER_no_eq, BER_raw_eq, BER_denoised_eq]);
set(gca, 'YScale', 'log');
b_ber.FaceColor = 'flat';
b_ber.CData = b_snr.CData;
% 1. Set the visual anchor point for the bottom of the bars
b_ber.BaseValue = 1e-8; % that way the bar is not drawn until -infinity
% 2. Manually set the Y-axis limits to add vertical padding
% Note: If your BER ever drops below 1e-4, just change this lower limit to 1e-5 or 1e-6
ylim([min(BER_raw_eq,min(BER_denoised_eq,BER_no_eq))/2, 1]);
%xlabel('Equalization Method');
ylabel('BER');
title('BER Comparison');
xticks(1:3);
xticklabels({'No EQ', 'Raw EQ', 'Denoised EQ'});
grid on;

while 1==1
% loop forever
end