Reliable cos^2(x) Approximation with Gaussian Peaks: Fix curve_fit by Dropping Peak Count and Seeding via find_peaks

Fit cos^2(x) with Gaussian peaks via scipy curve_fit: drop the non-differentiable peak count, estimate spacing and offsets with find_peaks to improve the fit.

Fitting a cos^2(x) waveform with a sum of Gaussian peaks looks straightforward, yet a naive application of scipy.optimize.curve_fit can lead to a painfully bad fit. The root cause is not the data or the function choice, but how the model parameters are defined and initialized.

Problem setup

The goal is to approximate cos^2(kx) with a stack of Gaussians spaced by a constant offset. The code below builds the target curve, a Gaussian basis, and attempts to fit all parameters at once, including how many peaks to use.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
freq = 1
mult = 3
x_lo = -mult * np.pi / (2 * freq)
x_hi =  mult * np.pi / (2 * freq)
grid = np.linspace(x_lo, x_hi, 1000)
target = np.cos(freq * grid) ** 2
def bell(xv, center, sigma):
    return np.exp(- (xv - center) ** 2 / (2 * sigma ** 2))
def bell_stack(xv, center, sigma, spacing, count):
    total = 0
    for idx in range(int(count)):
        pos = center + idx * spacing
        total += bell(xv, pos, sigma)
    return total
est = curve_fit(bell_stack, grid, target)[0]
c0 = est[0]
sg = est[1]
gap = est[2]
num = int(est[3])
approx = 0
for idx in range(num):
    pos = c0 + idx * gap
    approx += bell(grid, pos, sg)
fig, ax = plt.subplots()
ax.plot(grid, target, label="cosine function")
ax.plot(grid, approx, label="gaussian approximation to cosine")
ax.set_xlabel("x")
ax.set_ylabel("Amplitude")
plt.grid()
plt.legend()
plt.show()

Why the fit collapses

Two aspects work against the solver. First, the model includes a parameter that represents how many Gaussians to stack. That parameter is inherently discrete and non-differentiable. If an optimizer perturbs that value by an infinitesimal amount, nothing in the output changes. A local optimizer such as curve_fit relies on local gradient information, so it cannot learn a meaningful update for such a parameter. Second, the objective surface has many local minima. A local optimizer can get trapped early, and because it starts from the default initial guess of ones, this is especially harmful for the peak count and spacing.

There are also alternative modeling directions one might consider in other contexts, such as working in the Fourier domain or using the actual squared cosine instead of Gaussians, but here the requirement is to stay with Gaussian peaks in x-space.

Practical fix: remove the non-differentiable parameter and seed sensible starting values

The key is to take the discrete parameter out of the fit and to provide informed initial guesses. A reliable way to do that here is to detect the number and positions of peaks in the target signal. With scipy.signal.find_peaks, you can infer how many peaks to use, estimate the spacing between them, and pick the location of the first peak. Then, pass a model with a fixed peak count to curve_fit and let it only optimize continuous parameters.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
from functools import partial
freq = 1
mult = 3
x_lo = -mult * np.pi / (2 * freq)
x_hi =  mult * np.pi / (2 * freq)
grid = np.linspace(x_lo, x_hi, 1000)
target = np.cos(freq * grid) ** 2
def bell(xv, center, sigma):
    return np.exp(- (xv - center) ** 2 / (2 * sigma ** 2))
def bell_stack(xv, center, sigma, spacing, count):
    total = 0
    for idx in range(int(count)):
        pos = center + idx * spacing
        total += bell(xv, pos, sigma)
    return total
peak_idx, _ = find_peaks(target)
num = len(peak_idx)
init_gap = np.mean(np.diff(grid[peak_idx]))
init_c0 = grid[peak_idx[0]]
model = partial(bell_stack, count=num)
fit_params = curve_fit(model, grid, target, p0=(init_c0, 1, init_gap))[0]
c0, sg, gap = fit_params
approx = model(grid, *fit_params)
fig, ax = plt.subplots()
ax.plot(grid, target, label="cosine function")
ax.plot(grid, approx, label="gaussian approximation to cosine")
ax.set_xlabel("x")
ax.set_ylabel("Amplitude")
plt.grid()
plt.legend()
plt.show()

This approach fixes the non-differentiable degree of freedom, aligns the initial spacing with the detected waveform, and anchors the first peak. The optimizer is then free to refine continuous parameters rather than flailing on an integer.

Why this matters

Even with a reasonable model, optimization fails when parameterizations violate the assumptions of the solver. curve_fit is a local method; it needs gradients and decent seeds. Introducing discrete choices into a gradient-based routine or relying on default initial guesses invites poor convergence. In contrast, pre-estimating structural aspects of the signal and constraining the fit to continuous variables gives the algorithm something it can actually optimize.

Conclusion

When approximating structured signals like cos^2(x) with Gaussian sums, think carefully about what the optimizer can handle. Keep the parameter space differentiable by fixing discrete counts up front, and seed the fit with information extracted from the data, such as peak count, spacing, and a first-peak location. With those adjustments, the same fitting tools behave far more predictably and deliver the approximation you actually need.