The Quantum Fourier Transform: Build It Yourself in Qiskit

The Quantum Fourier Transform is the single most important subroutine in quantum computing. It's the engine inside Shor's factoring algorithm, quantum phase estimation, and order finding. It's also the most commonly misunderstood algorithm in the field — because on paper it looks like an exponentially faster FFT, and it absolutely is not.

This tutorial builds the QFT from scratch in Qiskit, proves the implementation is correct against Qiskit's own, and then explains carefully where the speedup goes.

All code in this post was executed against Qiskit 2.5.0 and Qiskit Aer 0.17.2 on Python 3.11. The outputs shown are real.

What the QFT actually does

The classical Discrete Fourier Transform takes a vector of N numbers and returns N Fourier coefficients. The QFT does exactly the same linear map — but the input vector is the amplitude vector of an n-qubit state, so N = 2ⁿ.

On basis states it acts as:

QFT |x⟩ = (1/√N) Σ_{k=0}^{N-1} e^{2πi·xk/N} |k⟩

Notice what's on the right: every basis state gets amplitude of the same magnitude (1/√N), and all the information about x lives in the complex phase e^{2πi·xk/N}. The QFT takes information that was sitting in the computational basis and pushes it into the phases.

That single fact is the whole story of the QFT. It's a perfectly uniform state as far as any measurement in the computational basis can tell. If you QFT a basis state and measure, you get a uniformly random n-bit string. Every time. You learn nothing.

The readout problem, stated plainly

Here's the trap. The QFT on n qubits takes O(n²) gates. Classically, an FFT on the same 2ⁿ amplitudes takes O(n·2ⁿ) operations. That is a genuine exponential gap in gate count, and it's the number everybody quotes.

But it isn't a usable speedup, for two reasons:

1. You can't load the input. Getting 2ⁿ arbitrary classical amplitudes into a quantum register generally costs O(2ⁿ) gates, wiping out the advantage before you start.
2. You can't read the output. The Fourier coefficients are amplitudes. You cannot look at amplitudes — you can only sample from their squared magnitudes, one n-bit string per shot. Extracting all 2ⁿ coefficients would take exponentially many repetitions.

So the QFT is never a drop-in replacement for the FFT. It is a transform you apply in the middle of a larger circuit, where the goal isn't to read the spectrum but to convert a periodicity that was hidden in phases into a peak you can actually sample. That's the shape every real application takes.

The circuit: Hadamard, rotate, recurse, swap

The construction follows directly from factoring the exponential above into a product over qubits. For the most significant qubit:

1. Apply a Hadamard. This creates a superposition whose relative phase encodes one bit of the input.
2. Apply controlled-phase gates CP(π/2^k) from each less significant qubit, controlled onto that same qubit. Qubit distance k contributes a rotation of π/2^k — the further away, the smaller the correction.
3. Recurse on the remaining n−1 qubits.
4. Swap the qubit order at the end.

Why the swaps? The recursion naturally produces the output register in reverse bit order. Qubit 0 ends up holding what should be the most significant output bit, and vice versa. This is a real bug, not a bookkeeping convention — if you feed an unswapped QFT into phase estimation you get a bit-reversed answer. The swaps fix it. (In hardware-aware compilation you sometimes do drop the swaps and relabel the wires downstream instead, which is free. But if your QFT is a standalone reusable block, keep them.)

Building it in Qiskit

Make sure you have Qiskit and the Aer simulator installed:

pip install qiskit qiskit-aer

Now the construction, written recursively so it mirrors the math:

import numpy as np
from qiskit import QuantumCircuit, transpile
from qiskit.quantum_info import Statevector, Operator
from qiskit_aer import AerSimulator


def qft_rotations(circuit, n):
    """Hadamard + controlled phases on the top qubit, then recurse."""
    if n == 0:
        return circuit
    n -= 1                      # index of the most significant qubit
    circuit.h(n)
    for qubit in range(n):
        # qubit is (n - qubit) positions below the target
        circuit.cp(np.pi / 2 ** (n - qubit), qubit, n)
    qft_rotations(circuit, n)   # recurse on the rest
    return circuit


def swap_registers(circuit, n):
    """Reverse the qubit order to undo the bit reversal."""
    for qubit in range(n // 2):
        circuit.swap(qubit, n - qubit - 1)
    return circuit


def qft(n):
    qc = QuantumCircuit(n, name="QFT")
    qft_rotations(qc, n)
    swap_registers(qc, n)
    return qc


my_qft = qft(4)
print(my_qft.draw(output="text"))
print(my_qft.count_ops())

Real output:

                                                                          ┌───┐»
q_0: ──────■───────────────────────────────■──────────────────────■───────┤ H ├»
           │                               │                ┌───┐ │P(π/2) └───┘»
q_1: ──────┼────────■──────────────────────┼────────■───────┤ H ├─■─────────X──»
           │        │                ┌───┐ │P(π/4)  │P(π/2) └───┘           │  »
q_2: ──────┼────────┼────────■───────┤ H ├─■────────■───────────────────────X──»
     ┌───┐ │P(π/8)  │P(π/4)  │P(π/2) └───┘                                     »
q_3: ┤ H ├─■────────■────────■─────────────────────────────────────────────────»
     └───┘

OrderedDict([('cp', 6), ('h', 4), ('swap', 2)])

Four Hadamards, six controlled-phase gates, two swaps. That matches the count: n Hadamards and n(n−1)/2 = 6 rotations. The angles descend π/2, π/4, π/8 as the control qubit gets further from the target — exactly as expected.

Proving it's correct

Never trust a hand-built QFT. Compare the full unitary against Qiskit's QFTGate:

from qiskit.circuit.library import QFTGate

n = 4
mine = qft(n)

reference = QuantumCircuit(n)
reference.append(QFTGate(n), range(n))

a, b = Operator(mine), Operator(reference)
print("Equivalent:", a.equiv(b))
print("Max element difference:", np.abs(a.data - b.data).max())

Real output:

Equivalent: True
Max element difference: 2.285583597498884e-15

Identical to floating-point precision — not merely equal up to a global phase. If you're on an older Qiskit, use from qiskit.circuit.library import QFT and QFT(n).decompose() instead.

Seeing it find a period

This is the one demonstration that shows what the QFT is for. Prepare a state that's a uniform superposition over multiples of 4 — a comb with period 4 — and transform it:

n = 4
psi = np.zeros(2 ** n, dtype=complex)
for x in range(0, 2 ** n, 4):   # |0> + |4> + |8> + |12>
    psi[x] = 1
psi /= np.linalg.norm(psi)

out = Statevector(psi).evolve(qft(n))
for i, p in enumerate(out.probabilities()):
    if p > 1e-6:
        print(f"|{i:0{n}b}> = {i:2d}   p = {p:.4f}")

Real output:

|0000> =  0   p = 0.2500
|0100> =  4   p = 0.2500
|1000> =  8   p = 0.2500
|1100> = 12   p = 0.2500

The input had period 4 in position space; the output is concentrated on multiples of 16/4 = 4 in frequency space. Sample this a few times, take the greatest common divisor of the outcomes, and you recover the period. That is the readable signal — not the coefficients, but a periodicity converted into a sampling distribution.

Confirm it on the simulator with actual shots:

qc = QuantumCircuit(n)
qc.initialize(psi, range(n))
qc.compose(qft(n), inplace=True)
qc.measure_all()

sim = AerSimulator()
counts = sim.run(transpile(qc, sim), shots=4096).result().get_counts()
print(sorted(counts.items(), key=lambda kv: -kv[1]))

Real output:

[('0100', 1079), ('1000', 1030), ('0000', 1018), ('1100', 969)]

Four peaks, roughly 1024 counts each. Nothing else appears.

Where the QFT earns its keep

Quantum phase estimation (QPE). Given a unitary U and an eigenstate |ψ⟩ with U|ψ⟩ = e^{2πiθ}|ψ⟩, QPE writes θ into a register's phases using controlled-U powers, then applies the inverse QFT to convert those phases into a binary number you can measure. Note the direction: QPE uses QFT⁻¹, because the phases are already there and the job is to pull them out. In Qiskit that's just qft(n).inverse().

Order finding. Given a and N, find the smallest r with a^r ≡ 1 (mod N). Modular exponentiation creates a register periodic in r; the QFT turns that period into measurable peaks, exactly as in the demo above.

Shor's algorithm. Factoring reduces to order finding, so Shor's is order finding plus classical post-processing (continued fractions and a GCD). The 1994 paper is worth reading precisely because the QFT is the only quantum part — everything else is number theory. This is also why the algorithm threatens RSA, and why post-quantum cryptography migration is happening now rather than later.

In all three cases the QFT is doing the same job: turning phase information into position information so a measurement can see it.

Approximate QFT: drop the tiny rotations

On real hardware, CP(π/2^k) for large k is a rotation smaller than your gate error. Executing it adds noise and buys nothing. The Approximate QFT (AQFT) simply skips every rotation with k above a cutoff, reducing the gate count from O(n²) to O(n log n):

def aqft(n, k_max):
    qc = QuantumCircuit(n, name="AQFT")
    for size in range(n, 0, -1):
        m = size - 1
        qc.h(m)
        for qubit in range(m):
            k = m - qubit
            if k <= k_max:
                qc.cp(np.pi / 2 ** k, qubit, m)
    for q in range(n // 2):
        qc.swap(q, n - q - 1)
    return qc

Measured fidelity against the exact 4-qubit QFT:

Keeping rotations down to π/4 recovers 98.6% fidelity with 17% fewer two-qubit gates. On larger registers the saving is dramatic: keeping only k ≤ O(log n) rotations is provably enough for Shor's algorithm to succeed, and on a noisy device the AQFT often outperforms the exact QFT because the errors you avoid outweigh the approximation you introduce. Always benchmark both on a simulator with a realistic noise model before committing.

Try this next

• Build qft(n).inverse() and use it to implement phase estimation on a T gate — you should recover θ = 1/8 exactly with 3 counting qubits.
• Run the 4-qubit QFT through transpile(qc, backend, optimization_level=3) for a real device and watch the CP gates decompose into CNOTs; the depth increase is sobering.
• Compare against a variational algorithm like QAOA to see the contrast between structured algorithms with proven speedups and heuristic NISQ ones.

The takeaway to keep: the QFT is fast, exact, and useless by itself. Its power is entirely in what you wrap around it.