An OFDM link with pilot-based channel estimation

What you will build

This project builds an OFDM link in one NumPy file and measures how well it copes with multipath. The transmitter places QPSK data on 48 subcarriers and known pilots on 16 more, runs a 64-point IFFT and adds a 16-sample cyclic prefix. The channel is an 8-tap Rayleigh fading filter, drawn afresh for every OFDM symbol, plus white Gaussian noise. The receiver removes the prefix, runs an FFT, estimates the channel at the pilots by least squares, interpolates to the data subcarriers and divides each subcarrier by its estimate.

The program compares four receivers: one that knows the channel exactly, and three that estimate it from pilots with different interpolators. It prints the mean squared error (MSE) of each channel estimate and the BER of each receiver from 0 dB to 30 dB SNR. In our run the receiver with perfect channel knowledge matched the textbook Rayleigh curve. Linear interpolation between pilots stalled at a BER of about 1.4 × 10⁻² because of interpolation error. DFT-based interpolation that keeps only the first 8 taps reached 7.42 × 10⁻⁴ at 30 dB, about 1.5 times the BER with perfect knowledge.

The OFDM from scratch and fading channels lessons cover the background. The QPSK modem project is the single-carrier version of the same measurement.

What you need

• Python 3.10 or later and NumPy. We ran it with Python 3.12 and NumPy 2.5.4.
• Nothing else. The whole run is vectorised and takes a few seconds.
• About 3 hours to read the design, run the program and try the variations.

Design

Numbers

The sample rate is 20 MHz and the FFT has 64 points, so the subcarrier spacing is 312.5 kHz. A symbol with its prefix is 80 samples, or 4.0 µs. The cyclic prefix is 16 samples (0.8 µs). Every 4th subcarrier, starting at subcarrier 0, carries a pilot of value 1, which gives 16 pilots. The other 48 subcarriers carry QPSK with unit average energy, the same as the pilots. Pilot overhead is 25 %.

Channel

The channel has 8 taps spaced one sample (50 ns) apart, so its maximum excess delay is 350 ns. The average tap powers fall off exponentially with a 1/e decay of 2 taps, normalised to a total of 1. That gives an rms delay spread of 82 ns. Each tap is a complex Gaussian with that average power, which makes the amplitude of each tap Rayleigh distributed. The program draws a new channel for every OFDM symbol. So each symbol must be estimated from its own pilots, with no averaging over time.

Because the channel is shorter than the prefix, the tail of each symbol falls inside the next symbol’s prefix, which the receiver discards. After the FFT, every subcarrier sees a single complex gain: Y[k] = H[k] X[k] + W[k]. That is the reason OFDM exists. A one-tap equaliser per subcarrier undoes the whole multipath channel.

SNR and the IFFT scale

The code uses a unitary IFFT and FFT, scaled by √64. Then the average power per time sample equals the average power per subcarrier, and white noise of variance N0 in time has variance N0 on every subcarrier. The SNR in the tables is Es/N0 per subcarrier. For QPSK, Eb/N0 is half of it, 3 dB lower.

Least-squares estimation

At pilot subcarrier k the receiver knows X[k] = 1, so the least-squares (LS) estimate is Y[k] / X[k]. It is unbiased, and its error is the noise on that subcarrier, with variance N0. LS on its own gives 16 noisy samples of H. The receiver needs H on all 64 subcarriers.

Interpolation

Linear. Draw straight lines between neighbouring pilots in the real and imaginary parts. Subcarriers 61 to 63 lie beyond the last pilot, so the code extends the last segment. Linear interpolation is cheap, but it cuts across the curvature of H. That leaves an error even with no noise.

DFT. The pilots sample H every 4 subcarriers. A 16-point IDFT of those samples returns the first 16 taps of the impulse response, exactly, provided the channel is no longer than 16 samples. Zero-padding those taps to 64 and running a 64-point FFT gives H on every subcarrier with no interpolation error at all. The cyclic prefix already limits useful channels to 16 samples, so this condition costs nothing.

DFT with truncation. If the receiver knows the channel is at most 8 taps long, it can zero taps 8 to 15 before the FFT. Those taps hold only noise, so the MSE halves. In practice the receiver would estimate the channel length or use the cyclic prefix length as a safe bound.

Theory for perfect channel knowledge

With perfect channel knowledge each data subcarrier is QPSK through a flat Rayleigh fading gain. The average BER per bit is 0.5 (1 − √(γ / (1 + γ))) with γ = Eb/N0. The program prints it next to the measured value as a check of the whole chain.

Code

Save this as ofdm_link.py. It is one block: the parameters, transmitter, channel, receiver, the three estimators and the sweep.

# ofdm_link.py: OFDM over a multipath channel with comb pilots, LS channel estimation,
# interpolation and one-tap equalisation. Prints channel-estimate MSE and BER against SNR.
import numpy as np

N_FFT = 64            # subcarriers
N_CP = 16             # cyclic prefix in samples
FS = 20e6             # sample rate, so subcarrier spacing is 312.5 kHz
PILOT_STEP = 4        # a pilot on every 4th subcarrier
N_TAPS = 8            # channel length in samples (0.4 us at 20 MHz)
DECAY_TAPS = 2.0      # exponential power delay profile, 1/e decay in samples
N_SYM = 20_000        # OFDM symbols per SNR point
rng = np.random.default_rng(11)

pilot_idx = np.arange(0, N_FFT, PILOT_STEP)                 # 16 pilots: 0, 4, ..., 60
data_idx = np.setdiff1d(np.arange(N_FFT), pilot_idx)        # 48 data subcarriers
pdp = np.exp(-np.arange(N_TAPS) / DECAY_TAPS)
pdp /= pdp.sum()                                            # unit average channel power


def qpsk(bits: np.ndarray) -> np.ndarray:
    b = bits.reshape(-1, 2)
    return ((1 - 2 * b[:, 0]) + 1j * (1 - 2 * b[:, 1])) / np.sqrt(2)


def random_channels(n: int) -> np.ndarray:
    """n independent Rayleigh channels, one per OFDM symbol, shape (n, N_TAPS)."""
    g = (rng.standard_normal((n, N_TAPS)) + 1j * rng.standard_normal((n, N_TAPS))) / np.sqrt(2)
    return g * np.sqrt(pdp)


def ofdm_tx(grid: np.ndarray) -> np.ndarray:
    """grid (n, N_FFT) of subcarrier values -> time samples with CP, shape (n, N_FFT + N_CP)."""
    x = np.fft.ifft(grid, axis=1) * np.sqrt(N_FFT)          # unitary IFFT: power per sample = per subcarrier
    return np.concatenate([x[:, -N_CP:], x], axis=1)


def multipath(x: np.ndarray, h: np.ndarray) -> np.ndarray:
    """Convolve each symbol with its own channel. The last N_TAPS - 1 output samples would land in
    the next symbol's cyclic prefix, which the receiver discards, so they are dropped here."""
    y = np.zeros((x.shape[0], x.shape[1] + N_TAPS - 1), dtype=complex)
    for k in range(N_TAPS):
        y[:, k:k + x.shape[1]] += h[:, [k]] * x
    return y[:, :x.shape[1]]


def ofdm_rx(y: np.ndarray) -> np.ndarray:
    return np.fft.fft(y[:, N_CP:], axis=1) / np.sqrt(N_FFT)


def interp_linear(h_p: np.ndarray) -> np.ndarray:
    """Linear interpolation between pilots, linear extrapolation past the last pilot."""
    k = np.arange(N_FFT)
    re = np.array([np.interp(k, pilot_idx, v) for v in h_p.real])
    im = np.array([np.interp(k, pilot_idx, v) for v in h_p.imag])
    out = re + 1j * im
    tail = k > pilot_idx[-1]                                # subcarriers 61, 62, 63
    slope = (h_p[:, -1] - h_p[:, -2]) / PILOT_STEP
    out[:, tail] = h_p[:, [-1]] + slope[:, None] * (k[tail] - pilot_idx[-1])
    return out


def interp_dft(h_p: np.ndarray, n_keep: int) -> np.ndarray:
    """DFT interpolation: pilots -> impulse response, keep n_keep taps, zero-pad, back to 64 bins."""
    g = np.fft.ifft(h_p, axis=1)                            # len(pilot_idx) taps
    g[:, n_keep:] = 0
    full = np.zeros((h_p.shape[0], N_FFT), dtype=complex)
    full[:, :len(pilot_idx)] = g
    return np.fft.fft(full, axis=1)


def run(snr_db: float) -> dict[str, tuple[float, float]]:
    n_bits = N_SYM * len(data_idx) * 2
    bits = rng.integers(0, 2, n_bits, dtype=np.int8)
    grid = np.zeros((N_SYM, N_FFT), dtype=complex)
    grid[:, data_idx] = qpsk(bits).reshape(N_SYM, -1)
    grid[:, pilot_idx] = 1.0                                # known pilot value, same power as data
    h = random_channels(N_SYM)
    H = np.fft.fft(h, n=N_FFT, axis=1)                      # true frequency response
    n0 = 10 ** (-snr_db / 10)                               # SNR = Es/N0 per subcarrier
    y = multipath(ofdm_tx(grid), h)
    y += np.sqrt(n0 / 2) * (rng.standard_normal(y.shape) + 1j * rng.standard_normal(y.shape))
    Y = ofdm_rx(y)
    h_ls = Y[:, pilot_idx] / grid[:, pilot_idx]             # least squares at the pilots
    estimates = {"perfect CSI": H,
                 "LS + linear": interp_linear(h_ls),
                 "LS + DFT, 16 taps": interp_dft(h_ls, 16),
                 "LS + DFT, 8 taps": interp_dft(h_ls, N_TAPS)}
    out = {}
    for name, H_hat in estimates.items():
        Z = Y[:, data_idx] / H_hat[:, data_idx]             # one-tap zero-forcing equaliser
        rx_bits = np.column_stack([(Z.real < 0).ravel(), (Z.imag < 0).ravel()]).astype(np.int8).ravel()
        ber = np.count_nonzero(rx_bits != bits) / n_bits
        mse = np.mean(np.abs(H_hat - H) ** 2)
        out[name] = (mse, ber)
    return out


print(f"{N_FFT} subcarriers, CP {N_CP}, {len(pilot_idx)} pilots, {len(data_idx)} data, "
      f"spacing {FS / N_FFT / 1e3:.1f} kHz, symbol {(N_FFT + N_CP) / FS * 1e6:.1f} us")
tau = np.arange(N_TAPS) / FS
rms_ds = np.sqrt(np.sum(pdp * tau ** 2) - np.sum(pdp * tau) ** 2)
print(f"channel: {N_TAPS} taps, rms delay spread {rms_ds * 1e9:.0f} ns, max excess delay {(N_TAPS - 1) / FS * 1e9:.0f} ns")
# Interpolation error alone: feed the true channel at the pilots, no noise.
h = random_channels(N_SYM)
H = np.fft.fft(h, n=N_FFT, axis=1)
for name, H_hat in [("linear", interp_linear(H[:, pilot_idx])), ("DFT, 8 taps", interp_dft(H[:, pilot_idx], N_TAPS))]:
    print(f"noise-free interpolation MSE, {name}: {np.mean(np.abs(H_hat - H) ** 2):.2e}")

names = ["perfect CSI", "LS + linear", "LS + DFT, 16 taps", "LS + DFT, 8 taps"]
results = {snr_db: run(snr_db) for snr_db in range(0, 31, 3)}
print("channel estimate MSE")
print("SNR (dB)  " + "".join(f"{n:>19}" for n in names[1:]))
for snr_db, r in results.items():
    print(f"{snr_db:>8}  " + "".join(f"{r[n][0]:>19.2e}" for n in names[1:]))
print(f"BER ({N_SYM * len(data_idx) * 2} bits per point)")
print("SNR (dB)  " + "".join(f"{n:>19}" for n in names) + "   theory (Rayleigh)")
for snr_db, r in results.items():
    g = 10 ** (snr_db / 10) / 2                              # Eb/N0 = SNR / 2 for QPSK
    theory = 0.5 * (1 - np.sqrt(g / (1 + g)))
    print(f"{snr_db:>8}  " + "".join(f"{r[n][1]:>19.2e}" for n in names) + f"   {theory:>17.2e}")

Every operation works on a two-dimensional array, one row per OFDM symbol. The only Python loops are over the 8 channel taps and over the SNR points.

Run

python -m pip install numpy
python ofdm_link.py

The sweep simulates 20,000 OFDM symbols, 1.92 million bits, at each of 11 SNR points. It takes about 3 s.

Test

This is the output from our run, with the seed fixed at 11.

64 subcarriers, CP 16, 16 pilots, 48 data, spacing 312.5 kHz, symbol 4.0 us
channel: 8 taps, rms delay spread 82 ns, max excess delay 350 ns
noise-free interpolation MSE, linear: 2.71e-02
noise-free interpolation MSE, DFT, 8 taps: 2.50e-32
channel estimate MSE
SNR (dB)          LS + linear  LS + DFT, 16 taps   LS + DFT, 8 taps
       0             8.10e-01           1.00e+00           5.01e-01
       3             4.20e-01           5.03e-01           2.51e-01
       6             2.23e-01           2.52e-01           1.26e-01
       9             1.25e-01           1.25e-01           6.28e-02
      12             7.68e-02           6.31e-02           3.15e-02
      15             5.20e-02           3.15e-02           1.58e-02
      18             3.97e-02           1.58e-02           7.91e-03
      21             3.38e-02           7.93e-03           3.97e-03
      24             3.03e-02           3.99e-03           1.99e-03
      27             2.87e-02           2.00e-03           9.99e-04
      30             2.83e-02           9.99e-04           4.99e-04
BER (1920000 bits per point)
SNR (dB)          perfect CSI        LS + linear  LS + DFT, 16 taps   LS + DFT, 8 taps   theory (Rayleigh)
       0             2.12e-01           2.96e-01           3.11e-01           2.77e-01            2.11e-01
       3             1.47e-01           2.19e-01           2.33e-01           1.99e-01            1.47e-01
       6             9.20e-02           1.49e-01           1.57e-01           1.29e-01            9.21e-02
       9             5.29e-02           9.40e-02           9.63e-02           7.62e-02            5.31e-02
      12             2.87e-02           5.79e-02           5.43e-02           4.20e-02            2.88e-02
      15             1.51e-02           3.72e-02           2.93e-02           2.25e-02            1.51e-02
      18             7.74e-03           2.55e-02           1.53e-02           1.15e-02            7.74e-03
      21             3.97e-03           1.96e-02           7.76e-03           5.92e-03            3.92e-03
      24             1.92e-03           1.61e-02           3.88e-03           2.88e-03            1.98e-03
      27             9.75e-04           1.46e-02           1.99e-03           1.46e-03            9.95e-04
      30             4.94e-04           1.38e-02           1.02e-03           7.42e-04            4.99e-04

Check these points in your output.

Interpolation alone. With the true channel at the pilots and no noise, linear interpolation leaves an MSE of 2.71 × 10⁻². DFT interpolation leaves 2.50 × 10⁻³², which is rounding error. That single line explains the floor in both charts.

Perfect channel knowledge. The measured BER follows the Rayleigh theory column within the statistical spread at every point, for example 9.75 × 10⁻⁴ against 9.95 × 10⁻⁴ at 27 dB. If it does not, the noise scaling or the IFFT normalisation is wrong. Note how slowly the curve falls: one decade of BER per 10 dB of SNR. On a fading channel the errors come from the subcarriers that sit in deep fades.

DFT estimators. The MSE of the 16-tap DFT estimator equals N0 (1.00 at 0 dB, 9.99 × 10⁻⁴ at 30 dB), which is the LS noise passed through unchanged. Keeping 8 taps halves it to 4.99 × 10⁻⁴ at 30 dB. The BER penalty against perfect knowledge stays roughly constant at high SNR: 7.42 × 10⁻⁴ against 4.94 × 10⁻⁴ at 30 dB for 8 taps, and 1.02 × 10⁻³ for 16 taps.

Linear interpolation. At low SNR linear interpolation averages noise between two pilots and beats the 16-tap DFT estimator (MSE 0.810 against 1.00 at 0 dB). Above about 12 dB its interpolation error dominates. The MSE levels off near 2.8 × 10⁻² and the BER near 1.4 × 10⁻², however much power you add.

On a real radio

Over the air, three things happen before channel estimation that this simulation skips.

• Timing. The receiver must find where each OFDM symbol starts. Any start point inside the part of the prefix that is free of the previous symbol’s tail works, which gives OFDM a tolerance that single-carrier links lack. A start point too late lets the next symbol leak in. The usual tool is a correlation against the repeated prefix or a training symbol with two identical halves.
• Carrier frequency offset. An offset of a fraction of the 312.5 kHz spacing breaks the orthogonality of the subcarriers and adds inter-carrier interference that no one-tap equaliser can remove. Estimate it from the phase rotation between repeated parts of a training symbol and correct it in time before the FFT. See the glossary entry for carrier frequency offset.
• A channel that moves. Here every symbol has its own independent channel. A real channel changes over its coherence time, so receivers average pilots over several symbols, and pilots are spread in both time and frequency. 3GPP TS 38.211 (NR physical channels and modulation), in the Release 19 version in our reference library, defines the demodulation reference signals that 5G NR uses for this job in clause 7.4.1.1.

Real systems also leave the DC subcarrier and band edges empty. DFT interpolation then works on an incomplete set of pilots, and truncation leaks energy at the edges. Expect the clean results above to degrade, and test edge subcarriers separately.

To try the link over the air, write the time samples as complex64, transmit them through an SDR at 20 MHz sample rate (or scale every number down to a rate your radio supports), and record the result. Use a cable with attenuators or a dummy load unless you hold a licence for the frequency. Add a training symbol at the start of each frame for the timing and frequency estimates. The lesson on recording and reading IQ covers the file side.

Things to try

1. Space the pilots out. Set PILOT_STEP = 8. Now 8 pilots sample H, and a channel of 8 taps sits exactly at the limit of DFT interpolation. In our run DFT interpolation was still exact (BER 1.05 × 10⁻³ at 30 dB), while the noise-free MSE of linear interpolation rose to 0.244 and its BER floor to about 9 × 10⁻². Then also set N_TAPS = 12. The impulse response now aliases in the 8-point IDFT, the noise-free MSE of the DFT estimator becomes 3.19 × 10⁻², and its BER floors near 1.8 × 10⁻².
2. Longer delay spread. Set DECAY_TAPS = 4.0 (rms delay spread 104 ns). H varies faster across frequency. In our run the noise-free MSE of linear interpolation rose from 2.71 × 10⁻² to 6.80 × 10⁻² and its BER at 30 dB was 2.93 × 10⁻². The DFT estimators did not care, because the channel still fits in 8 taps.
3. Break the prefix. Set N_CP = 4 with 8 channel taps. To see the effect, change multipath to process a continuous stream so the tail of one symbol really lands in the next. The BER floors even with perfect knowledge.
4. MMSE equalisation. Replace the zero-forcing division with H* Y / (|H|² + N0). For QPSK with hard decisions on a per-subcarrier basis the decisions do not change, because the scaling is positive. It starts to matter for 16-QAM, where the amplitude carries information.
5. 16-QAM. Map four bits per subcarrier and slice on a 4 × 4 grid. Now the amplitude of the channel estimate matters, and the gap between estimators widens.
6. Average over time. Hold each channel for 10 symbols and average the LS estimates over them. The MSE should drop by a factor of 10 for the DFT estimators.

Summary

• A cyclic prefix longer than the channel turns multipath into one complex gain per subcarrier. One division per subcarrier equalises it.
• LS estimation at the pilots costs nothing but leaves the full noise variance N0 in each estimate.
• Linear interpolation leaves an interpolation error (MSE 2.71 × 10⁻² here) that floors the BER near 1.4 × 10⁻² whatever the SNR.
• DFT interpolation is exact when the channel is shorter than the number of pilots. Truncating to the known channel length halved the estimation noise in this run.
• Measure the perfect-knowledge curve first. When it matches the Rayleigh formula, the transmitter, channel and noise scaling are right, and any remaining gap belongs to the estimator.

Background: the OFDM and multicarrier topic hub and the OFDM and NR numerology cheat sheet.