How to Zero-Pad a Spectrum Without Phase Shifts: Fix Discrete Fourier Series Indexing and DC Alignment

Learn why zero-padding the spectrum broke a rect signal: DFS indexing mismatch, DC bin misalignment, and wrong 2D vs 2D+1 scaling, plus how to get five copies.

Zero-padding a spectrum to stretch a discrete-time signal sounds straightforward, until index arithmetic bites. Here’s a compact walkthrough of why inserting zeros between Fourier series coefficients produced a complex, shifted waveform instead of five clean copies of a real rect signal—and how to fix it by getting the indexing and frequency scale right.

Problem setup

The task starts with a basic rect signal a[n] over the domain [-1000, 1000], equal to 1 for |n| < 100 and 0 elsewhere. Its discrete Fourier series coefficients are then computed. Next, the spectrum is “stretched” by inserting four zeros after every coefficient and scaling by 0.2, i.e., 0.2·[a0, 0, 0, 0, 0, a1, 0, 0, 0, 0, a2, …]. The expected result in time is a real signal made of five identical copies of the original, spaced apart.

Code that demonstrates the issue

The following snippet reproduces the behavior where the inverse transform returns a complex waveform and copies appear shifted. Names are illustrative, but the computational logic follows the described approach.

import numpy as np
import cmath
L = 1000
JJ = complex(0, 1)
PI = np.pi
M = 2 * L + 1
# Rect signal: 1 for |n| < 100
rect_sig = np.zeros(M)
for t in range(-99, 100):
    rect_sig[t + L] = 1
# Snap tiny FP noise and near-integers
eps_cut = 1e-10
def snap_near_integers(arr, tol=eps_cut):
    r = np.real(arr)
    im = np.imag(arr)
    r[np.abs(r) < tol] = 0
    im[np.abs(im) < tol] = 0
    r_fraction = r - np.round(r)
    im_fraction = im - np.round(im)
    r[np.abs(r_fraction) < tol] = np.round(r[np.abs(r_fraction) < tol])
    im[np.abs(im_fraction) < tol] = np.round(im[np.abs(im_fraction) < tol])
    return r + 1j * im
# Forward/inverse Discrete Fourier series (problematic scaling)
def fs_transform_bad(data, halfspan, inv=False):
    win_len = 2 * halfspan + 1  # used both as array length and frequency scale
    out = np.zeros(win_len, dtype=complex)
    if inv:
        for n in range(-halfspan, halfspan + 1):
            for k in range(-halfspan, halfspan + 1):
                out[n + halfspan] += data[k + halfspan] * cmath.exp(JJ * 2 * PI * k * n / win_len)
    else:
        for k in range(-halfspan, halfspan + 1):
            for n in range(-halfspan, halfspan + 1):
                out[k + halfspan] += (1 / win_len) * data[n + halfspan] * cmath.exp(-JJ * 2 * PI * k * n / win_len)
    return out
# Spectrum of rect
spec_a = fs_transform_bad(rect_sig, L)
spec_a = snap_near_integers(spec_a)
# Stretch in frequency by inserting zeros and scaling
factor = 5
spec_f = np.zeros(factor * M, dtype=complex)
spec_f[::factor] = 0.2 * spec_a
# Reconstruct with the same transform
guess_L = (len(spec_f) - 1) // 2
time_f = fs_transform_bad(spec_f, guess_L, inv=True)
time_f = snap_near_integers(time_f)

What went wrong and why it happens

The root cause is an indexing mismatch between D and N that breaks the position of the DC bin and shifts the spectral content. If the original half-range is D=L and the length is N=2D+1, multiplying N directly by a factor can make the new length even or, when odd, move the “middle” index away from DC. With factor 5 in this case, the array does have a middle, but the DC term no longer lands at the true center after the stride insertion. The zero-frequency coefficient that used to be at index L migrates by a couple of indices, which shifts the reconstructed signal and introduces an imaginary component.

You can see it directly by probing the padded spectrum. The center index of the stretched array is not where the DC sample of the original lands after zero insertion, so the inverse transform uses a misaligned DC bin. That misalignment alone explains the phase offset and the unexpected imaginary part even though the original a[n] is completely real and even.

There is a second, closely related issue in the discrete Fourier series implementation. The array length is 2D+1, but the frequency scale factor inside the exponent should be based on 2D, not 2D+1. Using 2D+1 both as array length and as the denominator in the angle moves the effective grid and compounds the shift when you later stretch. The size of the array remains 2D+1, but the normalization factor in the exponent must be 2D.

The fix: stretch D, not N, and use the correct frequency scale

Two targeted adjustments resolve the problem and produce five real, properly aligned copies of the signal. First, keep the invariant N = 2D + 1. When you stretch, multiply D by the factor, then reconstruct N from that D. Second, in the Fourier series code, use 2D as the normalization factor for the exponent while the array still has 2D+1 samples. After these corrections, the DC coefficient lands at the center index, and the inverse transform returns the expected real pattern.

import numpy as np
import cmath
L = 1000
JJ = complex(0, 1)
PI = np.pi
M = 2 * L + 1
# Rect signal
rect_sig = np.zeros(M)
for t in range(-99, 100):
    rect_sig[t + L] = 1
# FP cleanup
eps_cut = 1e-10
def snap_near_integers(arr, tol=eps_cut):
    r = np.real(arr)
    im = np.imag(arr)
    r[np.abs(r) < tol] = 0
    im[np.abs(im) < tol] = 0
    r_fraction = r - np.round(r)
    im_fraction = im - np.round(im)
    r[np.abs(r_fraction) < tol] = np.round(r[np.abs(r_fraction) < tol])
    im[np.abs(im_fraction) < tol] = np.round(im[np.abs(im_fraction) < tol])
    return r + 1j * im
# Correct Discrete Fourier series: array size is 2D+1, frequency scale is 2D
def fs_transform_ok(data, halfspan, inv=False):
    freq_span = 2 * halfspan            # scale for the exponent
    out = np.zeros(freq_span + 1, dtype=complex)  # array length still 2D+1
    if inv:
        for n in range(-halfspan, halfspan + 1):
            for k in range(-halfspan, halfspan + 1):
                out[n + halfspan] += data[k + halfspan] * cmath.exp(JJ * 2 * PI * k * n / freq_span)
    else:
        for k in range(-halfspan, halfspan + 1):
            for n in range(-halfspan, halfspan + 1):
                out[k + halfspan] += (1 / freq_span) * data[n + halfspan] * cmath.exp(-JJ * 2 * PI * k * n / freq_span)
    return out
# Forward transform
spec_a = fs_transform_ok(rect_sig, L)
spec_a = snap_near_integers(spec_a)
# Stretch by inserting zeros between coefficients, scaling by the factor
factor = 5
L_new = factor * L
M_new = 2 * L_new + 1
spec_g = np.zeros(M_new, dtype=complex)
spec_g[::factor] = spec_a / factor
# Inverse transform
time_g = fs_transform_ok(spec_g, L_new, inv=True)
time_g = snap_near_integers(time_g)

If you want to see the DC-bin misalignment from the broken version, a quick probe makes it obvious. The center entry is not the DC value after zero insertion; it sits a couple of indices away in this example with factor 5.

# Demonstrate the DC shift in the broken approach
factor = 5
spec_f = np.zeros(factor * M, dtype=complex)
spec_f[::factor] = 0.2 * spec_a
mid = spec_f.shape[0] // 2
print(spec_f[mid])       # ends up 0
print(spec_f[mid - 2])   # holds the scaled DC term

Why this matters

When manipulating Fourier series by index, the relationship between the half-range D, the length N = 2D + 1, and the DC bin is not optional bookkeeping; it defines the grid on which frequency operations occur. Stretching by touching N directly and mixing up the scale inside the exponential offsets the entire spectrum, so the inverse transform is forced to reconstruct from a shifted frequency grid. That explains the unexpected imaginary component and the loss of clean periodic copies in time.

Takeaways

Keep N = 2D + 1 as a strict invariant and stretch D, not N. In the discrete Fourier series, use 2D as the normalization factor in the exponent while the signal still has 2D + 1 samples. After inserting zeros between spectral samples, make sure the DC term lands at the true center index of the new spectrum. Do these three things and the inverse will yield the expected real, repeated waveform without spurious phase artifacts. If performance becomes a concern, consider vectorizing the inner sums, but correctness here is entirely about indexing and the frequency scale.