# Maximum of the sum of different sinusoids

How to find the maximum sum of several (more than two) harmonic sinusoidal oscillations of the form $$y_n(t) = A_n\sin(2\pi f_nt+\phi_n)$$ with different amplitudes $$A_n$$, frequencies $$f_n$$ and phase shifts $$\phi_n$$? Is there any formula or algorithm?

Not really.

If your frequencies are equidistant, i.e. something like $$f_n = n\cdot f_0$$, you can do an inverse FFT, which would be reasonably efficient.

Otherwise you can try to estimate it. The power of the signal will be

$$x_{rms} = \sqrt{\sum_n \frac{A_n^2}{2}} \tag{1}$$

If the number of frequencies is sufficiently large (> 10-ish) and the phases are sufficiently random, the amplitude probability distribution function will very much look like a Gaussian (Central Limit Theorem. While technically the maximum amplitude of a Gaussian signal is infinite, in practice you will find the peak to be between 3 and 5 times the RMS value.

To make things more complicated: this will depend on the length of the observation window. The longer the signal is observed, the higher the chances become that things stack up the wrong way (or right way)

The maximum possible amplitude is

$$x_{max} = \sum_n A_n \tag{2}$$

This happens when all individual frequencies line up in phase for a single time instance $$t_{max}$$ which would have to meet the condition (roughly speaking for cosines)

$$2\pi \cdot f_n \cdot t_{max} + \phi_n = 2 \pi M_n \forall n, M_n \in \mathbb{Z} \tag{3}$$

Since the peak of a cosine has a flat derivative we can soften this up to "close enough", i.e.

$$|2\pi \cdot f_n \cdot t_{max} + \phi_n - 2 \pi M_n| < \epsilon \forall n, M_n \in \mathbb{Z} \tag{4}$$

The smaller you make $$\epsilon$$ the closer you will get to the theoretical maximum but also the larger $$t_{max}$$ will become and it but it may a very large number indeed.

In any case you can bound the peak as

$$\sqrt{\sum_n \frac{A_n^2}{2}} \leq x_{peak} \leq \sum_n A_n \tag{4}$$

As Hilmar says, you can get pretty loose bounds.

Might be better just to grind it out:

import numpy as np
from matplotlib import pyplot as plt

def calculate_signal(A, omega, phi, t):
if isinstance(t, float):
signal = 0
else:
signal = np.zeros(len(t))
for index in range(len(A)):
signal = signal + A[index]*np.sin(omega[index]*t + phi[index])

return signal

A = [1,2,3,4]
omega = [ np.random.uniform(0, np.pi/4) for x in A]
phi = [np.random.uniform(0, np.pi*2) for x in A]

tt = np.linspace(0,100,10000)
xx = calculate_signal(A,omega,phi,tt)
plt.plot(tt,xx)

index_max = np.argmax(xx)

plt.plot(tt[index_max],xx[index_max],'r+')