Add DD-domain pilot pattern design for OTFS integrated comm/nav (LEO …

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

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

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)