diff --git a/applyChannelDD.m b/applyChannelDD.m new file mode 100644 index 0000000..f1a8772 --- /dev/null +++ b/applyChannelDD.m @@ -0,0 +1,63 @@ +function rxDD = applyChannelDD(txDD, chInfo, scs, cpSize, dopplerPrecomp_Hz) + +%-------------------------------------------------------------------------- +% +% Applies the multipath channel directly in the Delay-Doppler domain +% via 2D circular convolution. +% +% This is exact for any Doppler shift, unlike the TF-domain element-wise +% model (channelTF .* signal) which breaks down when fd/scs is large. +% +% Supports Doppler pre-compensation: subtracts a bulk Doppler shift +% (e.g. from satellite ephemeris) so the residual Doppler is small +% enough for the DD-domain pilot guard band to capture. +% +% The DD channel acts as: y[l,k] = sum_i h_i * x[(l-l_i) mod N, (k-k_i) mod M] +% where l_i and k_i are the delay and Doppler bin indices of path i. +% +%-------------------------------------------------------------------------- +% Input arguments: +% +% txDD N x M transmit DD grid +% chInfo Struct from multipathChannel with fields: +% .pathDelays_s - path delays in seconds +% .pathDopplers_Hz - path Doppler shifts in Hz +% .pathGains - path complex gains +% .numPaths - number of channel paths +% scs Subcarrier spacing (Hz) +% cpSize Cyclic prefix ratio +% dopplerPrecomp_Hz (optional) Bulk Doppler to pre-compensate (Hz) +% Typically from satellite ephemeris. Default = 0. +% +%-------------------------------------------------------------------------- +% Function returns: +% +% rxDD N x M received DD grid after channel (pre-compensated) +% +%-------------------------------------------------------------------------- + +if nargin < 5 + dopplerPrecomp_Hz = 0; +end + +[N, M] = size(txDD); + +% DD grid resolutions +Ts = (1 + cpSize) / scs; % OFDM symbol duration with CP (s) +delta_tau = 1 / (N * scs); % Delay resolution (s/bin) +delta_nu = 1 / (M * Ts); % Doppler resolution (Hz/bin) + +rxDD = zeros(N, M); + +for i = 1:chInfo.numPaths + % Map physical delay/Doppler to integer bin indices + % Subtract bulk Doppler pre-compensation from each path + li = round(chInfo.pathDelays_s(i) / delta_tau); % Delay bins + ki = round((chInfo.pathDopplers_Hz(i) - dopplerPrecomp_Hz) / delta_nu); % Residual Doppler bins + hi = chInfo.pathGains(i); + + % 2D circular convolution: shift and accumulate + rxDD = rxDD + hi * circshift(txDD, [li, ki]); +end + +end diff --git a/ddChannelEstimate.m b/ddChannelEstimate.m new file mode 100644 index 0000000..0bec727 --- /dev/null +++ b/ddChannelEstimate.m @@ -0,0 +1,181 @@ +function [hEst, delayEst, dopplerEst, navInfo] = ddChannelEstimate(rxGrid, pilotInfo) + +%-------------------------------------------------------------------------- +% +% Estimates the channel in the Delay-Doppler domain using the received +% pilot response. Also extracts navigation parameters (delay/Doppler) +% for LEO satellite integrated communication and navigation. +% +% Method: The transmitted pilot is a single impulse at (lp, kp). +% After passing through an L-path channel, the received DD grid shows +% the channel response as shifted/scaled copies of the pilot. +% The channel taps h_i at (l_i, k_i) can be read directly from the +% guard region of the received grid. +% +% Scanning is restricted to physically valid regions: +% - Delay: [0, maxDelayBins] (channel delays are non-negative) +% - Doppler: [-maxDopplerBins, +maxDopplerBins] (residual after pre-comp) +% This prevents data contamination from being detected as false paths. +% +%-------------------------------------------------------------------------- +% Input arguments: +% +% rxGrid N x M received DD grid (after SFFT demodulation) +% pilotInfo Struct from pilotPatternDD with pilot position info: +% .lp, .kp - pilot position +% .lGuard, .kGuard - guard band sizes (for data placement) +% .maxDelayBins - max physical delay (bins), scan bound +% .maxDopplerBins - max residual Doppler (bins), scan bound +% +%-------------------------------------------------------------------------- +% Function returns: +% +% hEst N x M estimated DD domain channel matrix +% delayEst Estimated delay indices of channel paths +% dopplerEst Estimated Doppler indices of channel paths +% navInfo Struct with navigation-relevant measurements +% +%-------------------------------------------------------------------------- + +[N, M] = size(rxGrid); +kp = pilotInfo.kp; +lp = pilotInfo.lp; +kGuard = pilotInfo.kGuard; +lGuard = pilotInfo.lGuard; + +% Physical scan bounds (restrict to valid channel response region) +if isfield(pilotInfo, 'maxDelayBins') + maxDelayBins = pilotInfo.maxDelayBins; +else + maxDelayBins = lGuard; % Fallback: scan full guard +end +if isfield(pilotInfo, 'maxDopplerBins') + maxDopplerBins = pilotInfo.maxDopplerBins; +else + maxDopplerBins = kGuard; % Fallback: scan full guard +end + +%% Define scan region (physically valid channel response area only) +% Delay: channel delays are non-negative, so scan [lp, lp + maxDelayBins] +% Doppler: after pre-comp, residual is bounded by ±maxDopplerBins +scanRows = lp : min(N, lp + maxDelayBins); +scanCols = max(1, kp - maxDopplerBins) : min(M, kp + maxDopplerBins); + +%% Detection threshold: pilot-relative +% Real channel responses have amplitude = pilot_amplitude * channel_gain. +% Use a threshold relative to the pilot to reject data contamination. +% With pilot boost of ~10 dB, pilot amplitude ≈ 3.16x data amplitude. +% Threshold at -8 dB below pilot catches paths with |h| > 0.4 (rel to LoS) +% while rejecting data contamination (amplitude ~ 1/pilot_boost). +pilotAmp = abs(rxGrid(lp, kp)); +threshRelativedB = -8; % Detect paths down to -8 dB below pilot +threshold = pilotAmp * 10^(threshRelativedB/20); + +%% Build DD channel impulse response +hEst = zeros(N, M); + +% Scan the valid region for significant taps +delayEst = []; +dopplerEst = []; +pathGains = []; + +for r = 1:length(scanRows) + for c = 1:length(scanCols) + lr = scanRows(r); + kc = scanCols(c); + val = rxGrid(lr, kc); + + if abs(val) > threshold + % Channel tap detected + delayShift = lr - lp; + dopplerShift = kc - kp; + + % Place in DD channel matrix (circular shift) + delayIdx = mod(delayShift, N) + 1; + dopplerIdx = mod(dopplerShift, M) + 1; + hEst(delayIdx, dopplerIdx) = val; + + delayEst = [delayEst; delayShift]; + dopplerEst = [dopplerEst; dopplerShift]; + pathGains = [pathGains; val]; + end + end +end + +%% If no paths detected, use the pilot position as single-tap channel +if isempty(delayEst) + hEst(1, 1) = rxGrid(lp, kp); + delayEst = 0; + dopplerEst = 0; + pathGains = rxGrid(lp, kp); +end + +%% Fractional Doppler estimation for LoS path (Quinn estimator) +% Quinn's second estimator uses complex DFT bins for optimal accuracy. +[~, losIdx] = max(abs(pathGains)); +losRow = lp + delayEst(losIdx); % row index in rxGrid +losCol = kp + dopplerEst(losIdx); % col index in rxGrid + +% Quinn estimator along Doppler (column) dimension +fracDopplerShift = dopplerEst(losIdx); % default: integer bin +if losCol > 1 && losCol < M + Xm1 = rxGrid(losRow, losCol - 1); + X0 = rxGrid(losRow, losCol); + Xp1 = rxGrid(losRow, losCol + 1); + + ap = real(Xp1 / X0); + am = real(Xm1 / X0); + dp = -ap / (1 - ap); + dm = am / (1 - am); + + if abs(dp) < abs(dm) + delta_k = dp; + else + delta_k = dm; + end + + delta_k = max(-0.5, min(0.5, delta_k)); + fracDopplerShift = dopplerEst(losIdx) + delta_k; +end + +% Quinn estimator along delay (row) dimension +fracDelayShift = delayEst(losIdx); % default: integer bin +if losRow > 1 && losRow < N + Xm1 = rxGrid(losRow - 1, losCol); + X0 = rxGrid(losRow, losCol); + Xp1 = rxGrid(losRow + 1, losCol); + + ap = real(Xp1 / X0); + am = real(Xm1 / X0); + dp = -ap / (1 - ap); + dm = am / (1 - am); + + if abs(dp) < abs(dm) + delta_l = dp; + else + delta_l = dm; + end + + delta_l = max(-0.5, min(0.5, delta_l)); + fracDelayShift = delayEst(losIdx) + delta_l; +end + +%% Navigation information extraction +pilotEnergy = abs(rxGrid(lp, kp))^2; +noiseEst = median(abs(rxGrid(:)).^2) * 0.5; % Noise from full grid + +navInfo.pathDelays = delayEst; +navInfo.pathDopplers = dopplerEst; +navInfo.pathGains = pathGains; +navInfo.pilotPower = pilotEnergy; +navInfo.noiseEst = noiseEst; +if noiseEst > 0 + navInfo.snrPilot = 10*log10(pilotEnergy / noiseEst); +else + navInfo.snrPilot = Inf; +end +navInfo.numPathsDetected = length(delayEst); +navInfo.losFracDoppler = fracDopplerShift; +navInfo.losFracDelay = fracDelayShift; + +end diff --git a/equaliser.m b/equaliser.m index 076353d..535e6d9 100644 --- a/equaliser.m +++ b/equaliser.m @@ -60,19 +60,23 @@ h(:,j) = interp1(sample(:,j),interpPoints); end -% Keep only 1st, 7th and 14th samples per carrier -h1 = [h(:,1)';h(:,7)';h(:,14)']'; +% Keep only 1st, middle and last samples per carrier +midSym = ceil(ofdmSym/2); +h1 = [h(:,1)'; h(:,midSym)'; h(:,ofdmSym)']'; % Linearly Interpolate samples between symbols -interpPoints1 = ( 1 : (1/6) : 2 ); % Calc. points to interpolate -interpPoints2 = ( (2+(1/7)) : (1/7) : 3 ); +seg1Len = midSym; % number of symbols in first segment (1 to midSym) +seg2Len = ofdmSym - midSym; % number of symbols in second segment (midSym to end) +interpPoints1 = linspace(1, 2, seg1Len); +interpPoints2 = linspace(2, 3, seg2Len + 1); +interpPoints2(1) = []; % avoid duplicating midSym for i = 1:size(txSig_matrix,1) - h17(i,:) = interp1(h1(i,:),interpPoints1); % Interpolate 1-7 - h714(i,:) = interp1(h1(i,:),interpPoints2); % interpolate 7-14 + h_seg1(i,:) = interp1(h1(i,:), interpPoints1); + h_seg2(i,:) = interp1(h1(i,:), interpPoints2); end -% Concatenate matracies -h = [h17';h714']'; +% Concatenate segments +h = [h_seg1'; h_seg2']'; % Set eq of CP carriers to same as 1st carrier for k = 1:(noSamples-numSC) diff --git a/main.m b/main.m index 6fc350c..01ff323 100644 --- a/main.m +++ b/main.m @@ -33,11 +33,12 @@ M = 16; % Modulation alphabet k = log2(M); % Bits/symbol cpSize = 0.07; % OFDM cyclic prefix size -scs = 15e3; % Subcarrier spacing Hz -Bw = 10e6; % System bandwidth Hz -ofdmSym = 14; % No. OFDM symbols / subframe +scs = 120e3; % Subcarrier spacing Hz (120 kHz for LEO-NTN) +Bw = 100e6; % System bandwidth Hz +ofdmSym = 128; % No. OFDM symbols / subframe (increased for LEO Doppler) EbNo = (-3:1:30)'; % Range of energy/bit to noise power ratio -velocity = 120; % Velocity of mobile rx relative tx km/hr +velocity = 26000; % Velocity of mobile rx relative LEO satellite km/hr +fc = 2e9; % Carrier frequency (Hz) S-band for LEO - must match multipathChannel.m codeRate = 2/4; % FEC code rate used maxIterations = 25; % Set maximum no. of iterations for LDPC decoder totalBits = 1e6; % The approx. total no. of bits simulated @@ -65,8 +66,37 @@ ldpcDecoder.MaximumIterationCount = maxIterations; % Assign decoder's maximum iterations noCodedbits = size(parityCheck_matrix,2); % Find the Codeword length +%-------------------------------------------------------------------------- +% LEO Satellite Pilot Pattern Configuration (DD domain) +%-------------------------------------------------------------------------- +% Guard bands sized to match LEO channel delay/Doppler spread +% Must match multipathChannel.m parameters (500 ns delay spread cap) +delta_tau_s = 1 / (numSC * scs); % Delay resolution (s) +maxDelayspread_s = min(0.5 * cpSize / scs, 500e-9); % LEO capped delay spread +v_ms = velocity * 1e3 / 3600; % Velocity in m/s +fd_hz = v_ms * fc / physconst('LightSpeed'); % Max Doppler shift (Hz) +Ts_sym = (1 + cpSize) / scs; % OFDM symbol duration with CP +delta_nu_hz = 1 / (ofdmSym * Ts_sym); % Doppler resolution (Hz) + +% Doppler pre-compensation from satellite ephemeris (standard in NR-NTN) +% Removes bulk Doppler, leaving only residual scatter spread +fd_precomp_hz = fd_hz; % Pre-comp value (from ephemeris) +residualScatter_hz = 0.05 * fd_hz; % Residual scatter after pre-comp +residualScatter_bins = ceil(residualScatter_hz / delta_nu_hz); + +maxDelayBins = ceil(maxDelayspread_s / delta_tau_s); % Physical max delay (bins) + +pilotConfig.kp = ceil(ofdmSym/2); % Pilot Doppler index (center) +pilotConfig.lp = ceil((Bw/scs - 12)/2); % Pilot delay index (center) +pilotConfig.lGuard = 2 * maxDelayBins + 1; % Delay guard: 2x for data protection +pilotConfig.kGuard = 2 * residualScatter_bins + 2; % Doppler guard: residual only (pre-comp) +pilotConfig.maxDelayBins = maxDelayBins; % Scan bound: physical delay range +pilotConfig.maxDopplerBins = residualScatter_bins + 1; % Scan bound: residual Doppler range +pilotConfig.pilotBoostdB = 10; % Pilot power boost (dB) + % Create Vectors for storing error data berOFDM = zeros(length(EbNo),3); berCOFDM = zeros(length(EbNo),3); berOTFS = zeros(length(EbNo),3); berCOTFS = zeros(length(EbNo),3); +berOTFS_pilot = zeros(length(EbNo),3); berCOTFS_pilot = zeros(length(EbNo),3); errorStats_coded = zeros(1,3); errorStats_uncoded = zeros(1,3); for repetition=1:repeats % Repeat simulation multiple times with a unqique channel for each repeat @@ -74,9 +104,9 @@ % Generate and Encode data [dataIn, dataBits_in, codedData_in, packetSize, numPackets, numCB] = dataGen(k,numDC,ofdmSym,totalBits,codeRate,ldpcEncoder); - % Generate Rayleigh Fading Channel Impulse Response - txSig_size = zeros((numSC+cpLen),ofdmSym); % Assign the size of the channel - rayChan = multipathChannel(cpSize, scs, txSig_size, velocity); % Create fading channel impulse response + % Generate TF-domain channel matrix H(subcarrier, symbol) + tfGridSize = zeros(numSC, ofdmSym); % TF grid dimensions + [channelTF, chInfo] = multipathChannel(cpSize, scs, tfGridSize, velocity); % TF-domain channel + DD params % QAM Modulator qamTx = qammod(dataIn,M,'InputType','bit','UnitAveragePower',true); % Apply QAM modulation @@ -95,45 +125,44 @@ % Calculate SNR snr = EbNo + 10*log10(codeRate*k) + 10*log10(numDC/((numSC))); - % Multicarrier Modulation - frameBuffer = guardbandTx; % Create a 'buffer' so subframes can be individually modulated - txframeBuffer = []; % Initilaise matrix + % Multicarrier Modulation (OFDM) + % Apply channel in TF domain, then modulate to time domain + frameBuffer = guardbandTx; + txframeBuffer_faded = []; % Faded time-domain signals + txframeBuffer_clean = []; % Clean time-domain signals for w = 1:packetSize - ofdmTx = modOFDM(frameBuffer(:,1:ofdmSym),numSC,cpLen,ofdmSym); % Apply OFDM modulation to a subframe of data - frameBuffer(:, 1:ofdmSym) = []; % Delete modulated data from frameBuffer - txframeBuffer = [txframeBuffer;ofdmTx]; % Add modulated subframe to transmit buffer + subframeTF = frameBuffer(:,1:ofdmSym); % TF-domain subframe + fadedTF = channelTF .* subframeTF; % Apply channel in TF domain + ofdmTx_faded = modOFDM(fadedTF, numSC, cpLen, ofdmSym); % Faded → time domain + ofdmTx_clean = modOFDM(subframeTF, numSC, cpLen, ofdmSym); % Clean → time domain + frameBuffer(:, 1:ofdmSym) = []; + txframeBuffer_faded = [txframeBuffer_faded; ofdmTx_faded]; + txframeBuffer_clean = [txframeBuffer_clean; ofdmTx_clean]; end - - + + % Loop through different values of EbNo for m = 1:length(EbNo) % Loop through the of packets to be transmitted for j = 1:numPackets rxframeBuffer = []; % Initialise matrix - + % Transmit each subframe individually for u = 1:packetSize - - % Remove next subframe from the transmit buffer - txSig = txframeBuffer( ((u-1)*numel(ofdmTx)+1) : u*numel(ofdmTx) ); - - % Apply Channel to input signal - fadedSig = zeros(size(txSig)); % Pre-allocate vector size - for i = 1:size(txSig,1) % Perform elementwise... - for j = 1:size(txSig,2) % ...matrix multiplication - fadedSig(i,j) = txSig(i,j).*rayChan(i,j); - end - end - + + % Extract faded and clean subframes + fadedSig = txframeBuffer_faded( ((u-1)*numel(ofdmTx_faded)+1) : u*numel(ofdmTx_faded) ); + txSig = txframeBuffer_clean( ((u-1)*numel(ofdmTx_clean)+1) : u*numel(ofdmTx_clean) ); + % AWGN Channel release(awgnChannel); powerDB = 10*log10(var(fadedSig)); % Calculate Tx signal power noiseVar = 10.^(0.1*(powerDB-snr(m))); % Calculate the noise variance rxSig = awgnChannel(fadedSig,noiseVar); % Pass the signal through a noisy channel - + % Equalisation eqSig = equaliser(rxSig,fadedSig,txSig,ofdmSym); - + % Demodulation rxSubframe = demodOFDM(eqSig,cpLen,ofdmSym); % Apply OFDM demodulation rxframeBuffer = [rxframeBuffer';rxSubframe']'; % Store demodulated subframe in rx buffer @@ -183,45 +212,43 @@ % Calculate SNR snr = EbNo + 10*log10(codeRate*k) + 10*log10(numDC/((numSC))) + 10*log10(sqrt(ofdmSym)); - % Multicarrier Modulation - frameBuffer = guardbandTx; % Create a 'buffer' so subframes can be individually modulated - txframeBuffer = []; % Initilaise matrix + % Multicarrier Modulation (OTFS) + % Apply channel in TF domain between ISFFT and OFDM mod + frameBuffer = guardbandTx; + txframeBuffer_faded = []; + txframeBuffer_clean = []; for w = 1:packetSize - otfsTx = ISFFT(frameBuffer(:,1:ofdmSym)); % Apply OTFS modulation to a subframe of data - ofdmTx = modOFDM(otfsTx,numSC,cpLen,ofdmSym); % Apply OFDM modulation - frameBuffer(:, 1:ofdmSym) = []; % Delete modulated data from frameBuffer - txframeBuffer = [txframeBuffer;ofdmTx]; % Add modulated subframe to transmit buffer + otfsTF = ISFFT(frameBuffer(:,1:ofdmSym)); % DD → TF domain + fadedTF = channelTF .* otfsTF; % Apply channel in TF domain + ofdmTx_faded = modOFDM(fadedTF, numSC, cpLen, ofdmSym); % Faded → time domain + ofdmTx_clean = modOFDM(otfsTF, numSC, cpLen, ofdmSym); % Clean → time domain + frameBuffer(:, 1:ofdmSym) = []; + txframeBuffer_faded = [txframeBuffer_faded; ofdmTx_faded]; + txframeBuffer_clean = [txframeBuffer_clean; ofdmTx_clean]; end - + % Loop through different values of EbNo for m = 1:length(EbNo) % Loop through the of packets to be transmitted for j = 1:numPackets rxframeBuffer = []; % Initialise matrix - + % Transmit each subframe individually for u = 1:packetSize - - % Remove next subframe from the transmit buffer - txSig = txframeBuffer( ((u-1)*numel(ofdmTx)+1) : u*numel(ofdmTx) ); - - % Apply Channel to input signal - fadedSig = zeros(size(txSig)); % Pre-allocate vector size - for i = 1:size(txSig,1) % Perform elementwise... - for j = 1:size(txSig,2) % ...matrix multiplication - fadedSig(i,j) = txSig(i,j).*rayChan(i,j); - end - end - + + % Extract faded and clean subframes + fadedSig = txframeBuffer_faded( ((u-1)*numel(ofdmTx_faded)+1) : u*numel(ofdmTx_faded) ); + txSig = txframeBuffer_clean( ((u-1)*numel(ofdmTx_clean)+1) : u*numel(ofdmTx_clean) ); + % AWGN Channel release(awgnChannel); powerDB = 10*log10(var(fadedSig)); % Calculate Tx signal power noiseVar = 10.^(0.1*(powerDB-snr(m))); % Calculate the noise variance rxSig = awgnChannel(fadedSig,noiseVar); % Pass the signal through a noisy channel - + % Equalisation eqSig = equaliser(rxSig,fadedSig,txSig,ofdmSym); - + % Demodulation otfsRx = demodOFDM(eqSig,cpLen,ofdmSym); % Apply OFDM demodulation rxSubframe = SFFT(otfsRx); % Apply OTFS demodulation @@ -262,15 +289,226 @@ berCOTFS(m,:) = errorStats_coded; % Save coded BER data errorStats_uncoded = errorRate(codedData_in,codedData_out,1); % Reset the error rate calculator errorStats_coded = errorRate1(dataBits_in,dataBits_out,1); % Reset the error rate calculator - + end - + + +%-------------------------------------------------------------------------- +% OTFS with DD-Domain Pilot Pattern (LEO Satellite) +%-------------------------------------------------------------------------- + + % Calculate SNR (same as OTFS) + snr = EbNo + 10*log10(codeRate*k) + 10*log10(numDC/((numSC))) + 10*log10(sqrt(ofdmSym)); + + % --- Transmitter: Place pilot + data on DD grid --- + % For each subframe, create DD grid with pilot pattern + frameBuffer_coded = guardbandTx; % Use same guard-banded QAM data + + % Store pilot info for all subframes (same pattern each subframe) + pilotConfig.lp = ceil(size(guardbandTx,1)/2); % Update to match actual grid rows + pilotConfig.kp = ceil(ofdmSym/2); + + txDDGrids = cell(packetSize, 1); + rxDDGrids_faded = cell(packetSize, 1); + ddGridSize = [size(guardbandTx,1), ofdmSym]; + + for w = 1:packetSize + subframe = frameBuffer_coded(:, 1:ofdmSym); + frameBuffer_coded(:, 1:ofdmSym) = []; + + % Extract data symbols from subframe (column-major order) + dataSyms = subframe(:); + + % Create DD grid with pilot pattern + [ddGrid, dataIdx, pilotIdx, guardIdx, pilotInfoTx] = ... + pilotPatternDD(dataSyms, size(subframe,1), ofdmSym, pilotConfig); + + % Apply channel in DD domain with Doppler pre-compensation + txDDGrids{w} = ddGrid; + rxDDGrids_faded{w} = applyChannelDD(ddGrid, chInfo, scs, cpSize, fd_precomp_hz); + end + + % --- Receiver: Add noise in DD domain and estimate/equalize --- + for m = 1:length(EbNo) + for j = 1:numPackets + rxframeBuffer = []; + + for u = 1:packetSize + rxDD_faded = rxDDGrids_faded{u}; + + % Add AWGN directly in DD domain + % (SFFT/ISFFT are unitary transforms, so AWGN statistics are preserved) + sigPower = 10*log10(mean(abs(rxDD_faded(:)).^2)); + noiseVar = 10^(0.1*(sigPower - snr(m))); + noise = sqrt(noiseVar/2) * (randn(size(rxDD_faded)) + 1j*randn(size(rxDD_faded))); + rxDD = rxDD_faded + noise; + + % DD-domain channel estimation using pilot + [hEst, delayEst, dopplerEst, navInfo] = ddChannelEstimate(rxDD, pilotInfoTx); + + % DD-domain equalization + [eqDD, eqDataSyms] = otfsEqualiserDD(rxDD, hEst, dataIdx, pilotInfoTx); + + % Reconstruct subframe from equalized data at data positions + rxSubframe = zeros(ddGridSize); + if length(eqDataSyms) >= length(dataIdx) + rxSubframe(dataIdx) = eqDataSyms(1:length(dataIdx)); + end + + rxframeBuffer = [rxframeBuffer'; rxSubframe']'; + end + + % Remove guard band nulls (same as original OTFS path) + parallelRx = rxframeBuffer; + parallelRx((numDC/2)+1:(numDC/2)+11, :) = []; + parallelRx(1:1, :) = []; + qamRx = reshape(parallelRx, [numel(parallelRx), 1]); + + % Uncoded demodulation + dataOut = qamdemod(qamRx, M, 'OutputType', 'bit', 'UnitAveragePower', true); + codedData_out = randdeintrlv(dataOut, 4831); + codedData_out(numel(codedData_in)+1:end) = []; + errorStats_uncoded = errorRate(codedData_in, codedData_out, 0); + + % Coded demodulation + powerDB = 10*log10(var(qamRx)); + noiseVar = 10.^(0.1*(powerDB-(EbNo(m) + 10*log10(codeRate*k) - 10*log10(sqrt(numDC))))); + dataOut = qamdemod(qamRx, M, 'OutputType', 'approxllr', 'UnitAveragePower', true, 'NoiseVariance', noiseVar); + codedData_out1 = randdeintrlv(dataOut, 4831); + codedData_out1(numel(codedData_in)+1:end) = []; + + dataBits_out = []; + dataOut_buffer = codedData_out1; + for q = 1:numCB + dataBits_out = [dataBits_out; ldpcDecoder(dataOut_buffer(1:noCodedbits))]; + dataOut_buffer(1:noCodedbits) = []; + end + dataBits_out = double(dataBits_out); + errorStats_coded = errorRate1(dataBits_in, dataBits_out, 0); + end + + berOTFS_pilot(m,:) = errorStats_uncoded; + berCOTFS_pilot(m,:) = errorStats_coded; + errorStats_uncoded = errorRate(codedData_in, codedData_out, 1); + errorStats_coded = errorRate1(dataBits_in, dataBits_out, 1); + end + + % Print navigation info from last subframe + fprintf('\n--- LEO Navigation Info (from DD-domain pilot) ---\n'); + fprintf('Detected %d channel paths\n', navInfo.numPathsDetected); + fprintf('Pilot SNR: %.1f dB\n', navInfo.snrPilot); + fprintf('Pilot overhead: %.1f%%\n', pilotInfoTx.overheadPercent); + fprintf('Effective pilot boost: %.1f dB\n', pilotInfoTx.effectivePilotBoostdB); + for pp = 1:length(navInfo.pathDelays) + fprintf(' Path %d: delay=%+d bins, Doppler=%+d bins, |gain|=%.3f\n', ... + pp, navInfo.pathDelays(pp), navInfo.pathDopplers(pp), abs(navInfo.pathGains(pp))); + end + fprintf('---------------------------------------------------\n'); + + %---------------------------------------------------------------------- + % Velocity & Range Estimation from DD-Domain Pilot Detection + %---------------------------------------------------------------------- + % Physical resolutions of the OTFS DD grid + c = physconst('LightSpeed'); % Speed of light (m/s) + Ts_sym = (1 + cpSize) / scs; % OFDM symbol duration with CP (s) + N_delay = size(rxDD, 1); % Number of delay bins (= numSC) + M_doppler = size(rxDD, 2); % Number of Doppler bins (= ofdmSym) + + delta_tau = 1 / (N_delay * scs); % Delay resolution (s/bin) + delta_nu = 1 / (M_doppler * Ts_sym); % Doppler resolution (Hz/bin) + + fprintf('\n=== Velocity & Range Estimation ===\n'); + fprintf('System parameters:\n'); + fprintf(' Carrier frequency fc = %.2f GHz\n', fc/1e9); + fprintf(' Subcarrier spacing df = %.0f kHz\n', scs/1e3); + fprintf(' Delay bins N = %d\n', N_delay); + fprintf(' Doppler bins M = %d\n', M_doppler); + fprintf(' Delay resolution d_tau = %.2f ns (%.2f m)\n', delta_tau*1e9, delta_tau*c); + fprintf(' Doppler resolution d_nu = %.2f Hz (%.2f m/s)\n', delta_nu, delta_nu*c/fc); + + % --- Identify LoS path (strongest gain) --- + [~, losIdx] = max(abs(navInfo.pathGains)); + los_delay_bin = navInfo.pathDelays(losIdx); + los_doppler_bin = navInfo.pathDopplers(losIdx); + + % Use fractional Doppler/delay from parabolic interpolation + los_frac_doppler = navInfo.losFracDoppler; % Fractional Doppler (bins) + los_frac_delay = navInfo.losFracDelay; % Fractional delay (bins) + + % Convert LoS path to physical values (using fractional estimates) + % Add back the pre-compensated bulk Doppler (from ephemeris) + % NOTE: SFFT convention maps positive Doppler to negative DD bins + % Physical Doppler = -DD_Doppler_shift * delta_nu + tau_los = los_frac_delay * delta_tau; % Propagation delay (s) + nu_residual = -los_frac_doppler * delta_nu; % Residual Doppler (Hz) — negate for SFFT sign + nu_los = nu_residual + fd_precomp_hz; % Total Doppler = residual + bulk (Hz) + range_los = tau_los * c; % One-way range (m) + v_est = nu_los * c / fc; % Estimated radial velocity (m/s) + v_est_kmh = v_est * 3.6; % Convert to km/hr + + % True velocity for comparison + v_true_kmh = velocity; % From simulation parameter (km/hr) + v_true_ms = velocity * 1e3 / 3600; % m/s + fd_true = v_true_ms * fc / c; % True max Doppler shift (Hz) + + fprintf('\n Doppler pre-compensation = %.2f Hz (from ephemeris)\n', fd_precomp_hz); + fprintf(' Residual Doppler guard = %d bins (kGuard)\n', pilotConfig.kGuard); + fprintf('\n--- LoS Path (strongest, Path %d) ---\n', losIdx); + fprintf(' Delay = %+d bins (frac: %+.2f) => tau = %.2f ns => range = %.2f m\n', ... + los_delay_bin, los_frac_delay, tau_los*1e9, range_los); + fprintf(' Doppler = %+d bins (frac: %+.2f) => nu_residual = %.2f Hz\n', ... + los_doppler_bin, los_frac_doppler, nu_residual); + fprintf(' Total Doppler = residual + precomp = %.2f + %.2f = %.2f Hz\n', ... + nu_residual, fd_precomp_hz, nu_los); + fprintf(' v_est = %.2f m/s (%.2f km/hr)\n', v_est, v_est_kmh); + fprintf('\n--- Comparison with True Velocity ---\n'); + fprintf(' True velocity = %.2f m/s (%.2f km/hr)\n', v_true_ms, v_true_kmh); + fprintf(' True max Doppler fd = %.2f Hz\n', fd_true); + fprintf(' Estimated velocity = %.2f m/s (%.2f km/hr)\n', v_est, v_est_kmh); + fprintf(' Velocity error = %.2f m/s (%.2f km/hr)\n', abs(v_est - v_true_ms), abs(v_est_kmh - v_true_kmh)); + fprintf(' Relative error = %.1f%%\n', abs(v_est - v_true_ms)/v_true_ms * 100); + + % --- All paths: physical delay and Doppler --- + fprintf('\n--- All Detected Paths (physical values) ---\n'); + fprintf(' %-6s %-12s %-12s %-14s %-14s %-10s\n', ... + 'Path', 'Delay(ns)', 'Range(m)', 'Doppler(Hz)', 'Velocity(m/s)', '|Gain|'); + for pp = 1:length(navInfo.pathDelays) + tau_pp = navInfo.pathDelays(pp) * delta_tau; + nu_pp = -navInfo.pathDopplers(pp) * delta_nu + fd_precomp_hz; % Negate for SFFT sign + add bulk Doppler + range_pp = tau_pp * c; + v_pp = nu_pp * c / fc; + fprintf(' %-6d %-12.2f %-12.2f %-14.2f %-14.2f %-10.3f\n', ... + pp, tau_pp*1e9, range_pp, nu_pp, v_pp, abs(navInfo.pathGains(pp))); + end + fprintf('==========================================\n'); + end %-------------------------------------------------------------------------- % Figures %-------------------------------------------------------------------------- -% Plot BER / EbNo curves -plotGraphs(berOFDM, berCOFDM, berOTFS, berCOTFS, M, numSC, EbNo); +% Plot BER / EbNo curves (including OTFS-Pilot results) +plotGraphs(berOFDM, berCOFDM, berOTFS, berCOTFS, M, numSC, EbNo, berOTFS_pilot, berCOTFS_pilot); + +% Plot DD-domain pilot pattern visualization +figure; +ddVis = zeros(size(guardbandTx,1), ofdmSym); +ddVis(dataIdx) = 0.3; % Data positions (gray) +ddVis(guardIdx) = 0; % Guard band (black) +ddVis(pilotIdx) = 1; % Pilot (bright) +imagesc(ddVis); +colormap(hot); +colorbar; +title('Delay-Doppler Grid: Pilot Pattern with Guard Bands'); +xlabel('Doppler Index (OFDM symbols)'); +ylabel('Delay Index (Subcarriers)'); +hold on; +% Draw guard band rectangle +rectangle('Position', [pilotConfig.kp-pilotConfig.kGuard-0.5, ... + pilotConfig.lp-pilotConfig.lGuard-0.5, ... + 2*pilotConfig.kGuard+1, 2*pilotConfig.lGuard+1], ... + 'EdgeColor', 'c', 'LineWidth', 2, 'LineStyle', '--'); +legend('Guard Band Boundary'); +hold off; diff --git a/multipathChannel.m b/multipathChannel.m index c3dd9ad..94e4d7e 100644 --- a/multipathChannel.m +++ b/multipathChannel.m @@ -1,22 +1,25 @@ -function [outSig] = multipathChannel(cpSize, delta_f, inSig, velocity) +function [H, chInfo] = multipathChannel(cpSize, delta_f, inSig, velocity) %-------------------------------------------------------------------------- % -% Generates and encodes random binary data +% Generates the time-frequency (TF) domain channel transfer function +% for an LEO satellite Rician fading channel with Doppler shift. +% +% Returns a 2D matrix H(subcarrier, symbol) for element-wise +% multiplication with the TF-domain signal (after ISFFT, before OFDM mod). % %-------------------------------------------------------------------------- -% Input arguments: +% Input arguments: % -% cpSize Size of the cyclic prefix +% cpSize OFDM cyclic prefix ratio % delta_f OFDM subcarrier spacing (Hz) -% inSig Input signal dimensions used to generate channel dimensions -% totalBits The approximate total bits to be simulated by the system +% inSig Matrix whose size defines the TF grid (N_sc x M_sym) % velocity Channel mobility in km/hr -% +% %-------------------------------------------------------------------------- -% Function returns: -% -% outSig Channel impulse response (serial vector) +% Function returns: +% +% H N_sc x M_sym TF-domain channel matrix % %-------------------------------------------------------------------------- % @@ -29,74 +32,74 @@ % %-------------------------------------------------------------------------- +% Get TF grid dimensions +[N_sc, M_sym] = size(inSig); % subcarriers x OFDM symbols -% Create N x M channel matrix -[N, M]= size(inSig'); % Size of inSig is used to create channel model -n = zeros(1,N); % delay_Doppler rows (doppler) -m = zeros(1,M); % delay_Doppler cols (delay) -H = transpose(n).*m; % Create matrix - -% Generate Channel Parameter -maxDelayspread = 0.5*((cpSize)/delta_f); % Calculate max delay spread with 1/2 rule pf thumb -L = round(2*maxDelayspread * M*delta_f); % Calculate necesary paths from bandwidth +% Generate Channel Parameters (LEO satellite Rician channel) +maxDelayspread = min(0.5*((cpSize)/delta_f), 500e-9); % Cap at 500 ns for LEO +L = min(round(2*maxDelayspread * N_sc*delta_f), 6); % Cap at 6 paths for LEO +L = max(L, 2); % At least 2 paths step = maxDelayspread/L; % calculate difference between delays pathDelays = (0:step:maxDelayspread); % Discrete even delays of L-path channel range shuffle; % Shuffle random no. generator -avgPathGains_dB = -(randi([3,7],[L,1])); % Generate random path gains in dB -avgPathGains = 10.^(0.1*avgPathGains_dB); % Convert to linear + +% Rician K-factor for LEO satellite (strong LoS component) +K_ric_dB = 10; % K-factor in dB +K_ric = 10^(K_ric_dB/10); % Linear K-factor +avgPathGains = zeros(L, 1); +avgPathGains(1) = sqrt(K_ric / (K_ric + 1)); % LoS path (dominant) +% NLOS paths: exponentially decaying power profile +nlosPowerTotal = 1 / (K_ric + 1); +nlosDecay = exp(-(0:L-2) / max(1, (L-2)/2)); +nlosDecay = nlosDecay / sum(nlosDecay); % Normalize to sum=1 +for i = 2:L + avgPathGains(i) = sqrt(nlosPowerTotal * nlosDecay(i-1)); +end % Calculate Max Doppler Shift v = velocity*1e3/3600; % Mobile speed (m/s) -fc = 3.5e9; % Carrier frequency +fc = 2e9; % Carrier frequency (S-band for LEO) fd = round(v*fc/physconst('lightspeed'));% Maximum Doppler shift to nearest Hz -% Generate doppler spreads w/ Jake's model +% Generate Doppler shifts for LEO satellite channel +% All paths share the bulk Doppler shift fd from satellite motion. +% Multipath scattering adds a small random perturbation on top. +scatterSpread = 0.05 * fd; % scattering spread << fd (5% of fd) +Vi = zeros(1, L); for l=0:L-1 - Vi(l+1) = fd* cos( (2*pi*l)/(L-1) ); + Vi(l+1) = fd + scatterSpread * cos( (2*pi*l)/(L-1) ); end % Initialize channel variables T = 1/delta_f; % unextended OFDM symbol period -Ts = (1+cpSize)/delta_f; % OFDM symbol period) +Ts = (1+cpSize)/delta_f; % OFDM symbol period with CP Ti = pathDelays; % Path delays hi = avgPathGains; % Path gains -% Create matrix representation of channel -for m=1:M % Loop along the rows - for n=1:N % Loop down the cols - - for x=1:L % Loop to sum terms in the channel memory - % Define terms of model - expTerm = (-2*1i*(pi)) * ((m+M/2).*delta_f.*Ti(x) - Vi(x).*(n).*Ts); +% Return channel parameters for DD-domain processing +chInfo.pathDelays_s = Ti; % Path delays in seconds +chInfo.pathDopplers_Hz = Vi; % Path Doppler shifts in Hz +chInfo.pathGains = hi; % Path complex gains +chInfo.numPaths = L; % Number of paths +chInfo.fd = fd; % Max Doppler shift (Hz) + +% Create TF-domain channel matrix H(subcarrier, symbol) +% H(k, n) = sum_i hi(i) * (1 + j*pi*Vi(i)*T) * exp(-j2pi*(k_freq*df*Ti(i) - Vi(i)*n*Ts)) +% where k_freq is the centered subcarrier index +H = zeros(N_sc, M_sym); + +for sc = 1:N_sc % Loop over subcarriers (rows) + for sym = 1:M_sym % Loop over OFDM symbols (columns) + for x = 1:L % Sum over L multipath components + % Centered subcarrier frequency index + freqIdx = sc + N_sc/2; + % TF channel: delay causes frequency-dependent phase, + % Doppler causes time-dependent phase + expTerm = (-2*1i*(pi)) * (freqIdx.*delta_f.*Ti(x) - Vi(x).*(sym).*Ts); hiPrime = hi(x)*(1 + 1i*(pi).*Vi(x).*T); - % Produce channel impulse response model - H(n, m) = H(n, m) + exp(expTerm) * hiPrime; - + H(sc, sym) = H(sc, sym) + exp(expTerm) * hiPrime; end end - end - -% % 3D plot of channnel model -% figure -% absH = 10*log10(abs(H)); -% mesh(absH) -% surf(absH) -% colormap(jet) % change color map -% shading interp % interpolate colors across lines and faces -% hold on; -% title(['Time-Frequency Plot of The Fading Channel']); -% xlabel('Frequency (MHz)'); -% xticks([8 142 274 408 540]) -% xticklabels({'3.496','3.498','3.500','3.502','3.504'}) -% ylabel('Time (ms)'); -% yticks([0 14 28 42 56 70 84 98 112 126 140]) -% yticklabels({'0','1','2','3','4','5','6','7','8','9','10'}) -% zlabel('Power (dB'); -% hold off; - - -% Convert to serial -outSig = reshape(H',[n*m,1]); end diff --git a/otfsEqualiserDD.m b/otfsEqualiserDD.m new file mode 100644 index 0000000..4163c9e --- /dev/null +++ b/otfsEqualiserDD.m @@ -0,0 +1,72 @@ +function [eqGrid, dataSymbols] = otfsEqualiserDD(rxGrid, hEst, dataIdx, pilotInfo) + +%-------------------------------------------------------------------------- +% +% Equalizes the received OTFS signal in the Delay-Doppler domain +% using the estimated DD channel response. +% +% Method: Constructs the effective DD channel matrix and performs +% MMSE or ZF equalization to recover data symbols. +% +%-------------------------------------------------------------------------- +% Input arguments: +% +% rxGrid N x M received DD grid (after SFFT demodulation) +% hEst N x M estimated DD channel matrix (from ddChannelEstimate) +% dataIdx Linear indices of data positions (from pilotPatternDD) +% pilotInfo Struct with pilot metadata +% +%-------------------------------------------------------------------------- +% Function returns: +% +% eqGrid N x M equalized DD grid +% dataSymbols Column vector of equalized data symbols (at dataIdx) +% +%-------------------------------------------------------------------------- + +[N, M] = size(rxGrid); + +%% Method: Frequency-domain MMSE equalization +% Convert DD channel to TF domain, equalize there, convert back. +% This is more computationally efficient than building the full +% N*M x N*M DD domain channel matrix. + +% Transform received DD grid to TF domain +rxTF = sqrt(N/M) * fft(ifft(rxGrid, [], 1), [], 2); % ISFFT + +% Transform DD channel to TF domain +hTF = sqrt(N/M) * fft(ifft(hEst, [], 1), [], 2); % ISFFT + +% Estimate noise variance from guard region +kp = pilotInfo.kp; +lp = pilotInfo.lp; +kGuard = pilotInfo.kGuard; +lGuard = pilotInfo.lGuard; +guardRows = max(1, lp - lGuard) : min(N, lp + lGuard); +guardCols = max(1, kp - kGuard) : min(M, kp + kGuard); + +% Noise from guard positions that are not near detected paths +guardVals = rxGrid(guardRows, guardCols); +% Use the weakest values in guard region as noise estimate +sortedGuard = sort(abs(guardVals(:)).^2); +numNoisesamples = max(1, floor(length(sortedGuard) * 0.3)); +noiseVar = mean(sortedGuard(1:numNoisesamples)); +if noiseVar == 0 + noiseVar = 1e-10; % Prevent division by zero +end + +%% Per-subcarrier MMSE equalization in TF domain +eqTF = zeros(N, M); +for col = 1:M + H_col = hTF(:, col); % Channel frequency response for this symbol + % MMSE equalizer: W = H* / (|H|^2 + sigma^2) + eqTF(:, col) = conj(H_col) ./ (abs(H_col).^2 + noiseVar) .* rxTF(:, col); +end + +%% Convert equalized TF grid back to DD domain via SFFT +eqGrid = sqrt(M/N) * ifft(fft(eqTF, [], 1), [], 2); % SFFT + +%% Extract data symbols from equalized grid +dataSymbols = eqGrid(dataIdx); + +end diff --git a/pilotPatternDD.m b/pilotPatternDD.m new file mode 100644 index 0000000..e542604 --- /dev/null +++ b/pilotPatternDD.m @@ -0,0 +1,132 @@ +function [ddGrid, dataIdx, pilotIdx, guardIdx, pilotInfo] = pilotPatternDD(dataSymbols, N, M, pilotConfig) + +%-------------------------------------------------------------------------- +% +% Places pilot, guard bands, and data symbols on the Delay-Doppler grid +% for OTFS integrated communication and navigation (LEO satellite). +% +% Uses the "Impulse Pilot + Guard Band" scheme in the DD domain. +% The pilot is a single high-energy impulse at (kp, lp), surrounded +% by a rectangular zero-value guard region. Data symbols fill the rest. +% +%-------------------------------------------------------------------------- +% Input arguments: +% +% dataSymbols Column vector of QAM data symbols to place on the grid +% N Number of rows (delay bins / subcarriers) +% M Number of columns (Doppler bins / OFDM symbols) +% pilotConfig Struct with fields: +% .kp Pilot Doppler index (1-based), default = ceil(M/2) +% .lp Pilot delay index (1-based), default = ceil(N/2) +% .kGuard Doppler guard band width (one side), default = 3 +% .lGuard Delay guard band width (one side), default = 2 +% .pilotBoostdB Pilot power boost in dB relative to avg data power +% .totalPowerdBm (optional) Total transmit power constraint in dBm +% +%-------------------------------------------------------------------------- +% Function returns: +% +% ddGrid N x M matrix: the delay-Doppler grid with pilot+data +% dataIdx Linear indices of data symbol positions +% pilotIdx Linear index of the pilot position +% guardIdx Linear indices of guard band (zero) positions +% pilotInfo Struct with pilot/guard metadata for receiver +% +%-------------------------------------------------------------------------- + +%% Default pilot configuration +if nargin < 4 || isempty(pilotConfig) + pilotConfig = struct(); +end +if ~isfield(pilotConfig, 'kp'), pilotConfig.kp = ceil(M/2); end +if ~isfield(pilotConfig, 'lp'), pilotConfig.lp = ceil(N/2); end +if ~isfield(pilotConfig, 'kGuard'), pilotConfig.kGuard = 3; end +if ~isfield(pilotConfig, 'lGuard'), pilotConfig.lGuard = 2; end +if ~isfield(pilotConfig, 'pilotBoostdB'), pilotConfig.pilotBoostdB = 10; end + +kp = pilotConfig.kp; +lp = pilotConfig.lp; +kGuard = pilotConfig.kGuard; +lGuard = pilotConfig.lGuard; + +%% Validate guard band fits within the grid +if (kp - kGuard < 1) || (kp + kGuard > M) + error('Doppler guard band exceeds grid boundary. Reduce kGuard or adjust kp.'); +end +if (lp - lGuard < 1) || (lp + lGuard > N) + error('Delay guard band exceeds grid boundary. Reduce lGuard or adjust lp.'); +end + +%% Initialize grid +ddGrid = zeros(N, M); + +%% Identify guard region indices (rectangular block around pilot) +guardRows = (lp - lGuard):(lp + lGuard); % delay indices +guardCols = (kp - kGuard):(kp + kGuard); % Doppler indices + +% Create logical mask for guard+pilot region +guardMask = false(N, M); +guardMask(guardRows, guardCols) = true; + +% Pilot position +pilotMask = false(N, M); +pilotMask(lp, kp) = true; + +% Guard = guard region minus the pilot itself +guardOnlyMask = guardMask & ~pilotMask; + +% Data positions = everything outside the guard region +dataMask = ~guardMask; + +%% Get linear indices +pilotIdx = find(pilotMask); +guardIdx = find(guardOnlyMask); +dataIdx = find(dataMask); + +%% Calculate number of available data positions +numDataPositions = length(dataIdx); + +%% Place data symbols +% Truncate or zero-pad data to fit available positions +if length(dataSymbols) > numDataPositions + % More data than available slots: truncate + ddGrid(dataIdx) = dataSymbols(1:numDataPositions); +elseif length(dataSymbols) < numDataPositions + % Fewer data symbols: place what we have, rest stays zero + ddGrid(dataIdx(1:length(dataSymbols))) = dataSymbols; +else + ddGrid(dataIdx) = dataSymbols; +end + +%% Place pilot symbol with power boost +% Pilot amplitude = sqrt(boostLinear) * avg_data_amplitude +% This ensures pilot energy is boosted relative to data +avgDataPower = mean(abs(dataSymbols).^2); +pilotBoostLinear = 10^(pilotConfig.pilotBoostdB / 10); +pilotAmplitude = sqrt(pilotBoostLinear * avgDataPower); +ddGrid(lp, kp) = pilotAmplitude; % Real-valued impulse pilot + +%% Power normalization under total power constraint +% Total power = pilot power + data power +% Normalize so that total average power per symbol = 1 +totalPower = sum(abs(ddGrid(:)).^2); +numActiveSymbols = numDataPositions + 1; % data + pilot +if totalPower > 0 + normFactor = sqrt(numActiveSymbols / totalPower); + ddGrid = ddGrid * normFactor; +end + +%% Store pilot info for receiver side +pilotInfo.kp = kp; +pilotInfo.lp = lp; +pilotInfo.kGuard = kGuard; +pilotInfo.lGuard = lGuard; +pilotInfo.numDataPositions = numDataPositions; +pilotInfo.numGuardPositions = length(guardIdx); +pilotInfo.pilotBoostdB = pilotConfig.pilotBoostdB; +pilotInfo.overheadPercent = 100 * (length(guardIdx) + 1) / (N * M); +pilotInfo.pilotPower = abs(ddGrid(lp, kp))^2; +pilotInfo.avgDataPower = mean(abs(ddGrid(dataIdx)).^2); +pilotInfo.effectivePilotBoostdB = 10*log10(pilotInfo.pilotPower / pilotInfo.avgDataPower); + +end diff --git a/plotGraphs.m b/plotGraphs.m index f69a76f..a6e4c60 100644 --- a/plotGraphs.m +++ b/plotGraphs.m @@ -1,22 +1,25 @@ -function [] = plotGraphs(berOFDM, berCOFDM, berOTFS, berCOTFS, M, numSC, EbNo) +function [] = plotGraphs(berOFDM, berCOFDM, berOTFS, berCOTFS, M, numSC, EbNo, berOTFS_pilot, berCOTFS_pilot) %-------------------------------------------------------------------------- % % Plots and formats BER graphs of the input BER vectors % %-------------------------------------------------------------------------- -% Input arguments: -% +% Input arguments: +% % berOFDM OFDM ber vector % berCOFDM coded-OFDM ber vector % berOTFS OTFS ber vector -% berCOTFS coded-OFDM ber vector +% berCOTFS coded-OTFS ber vector % M Modulation order % numSC no. subcarriers used in system +% EbNo Eb/N0 range +% berOTFS_pilot (optional) OTFS with DD pilot ber vector +% berCOTFS_pilot (optional) coded-OTFS with DD pilot ber vector % %-------------------------------------------------------------------------- -% Function returns: -% +% Function returns: +% % void % %-------------------------------------------------------------------------- @@ -34,7 +37,7 @@ berCOFDM(~berCOFDM)=1e-12; berCOTFS(~berCOTFS)=1e-12; -% Plotting +% Plotting figure % Plot non-coded results @@ -44,13 +47,22 @@ % Plot coded results semilogy(EbNo,berCOFDM(:,1),'--g'); %Plot simulated BER w/ C-OFDM semilogy(EbNo,berCOTFS(:,1),'--r'); %Plot simulated BER w/ C-OTFS - + +% Plot OTFS-Pilot results if provided +legendEntries = {'OFDM', 'OTFS', 'C-OFDM', 'C-OTFS'}; +if nargin >= 9 && ~isempty(berOTFS_pilot) + berCOTFS_pilot(~berCOTFS_pilot)=1e-12; + semilogy(EbNo,berOTFS_pilot(:,1),'-b','LineWidth',1.5); + semilogy(EbNo,berCOTFS_pilot(:,1),'--b','LineWidth',1.5); + legendEntries = [legendEntries, {'OTFS-Pilot (DD)', 'C-OTFS-Pilot (DD)'}]; +end + % Formatting graph axis([0 30 0.0005 0.5]) title(['BER for ', num2str(M), '-QAM in Wideband Rayleigh Fading with ', num2str(numSC),' Subcarriers']); ylabel('Bit Error Rate'); xlabel('EbNo (dB)'); -legend('OFDM', 'OTFS','C-OFDM','C-OTFS'); +legend(legendEntries); grid on; hold off;