# My summed sin waves seems to only work if my wavelength is divisible with the time steps. A problem with phase likely

I have created a program in python that can create any generic signal and a time range to display it as shown by the red '--' line below. I then created a way to calculate the phase, amplitude, and wavelength for each frequency in this signal via fft. Then, I create an equation for each frequency based on these amp, ph, and wavelengths. I then sum all of these equations and the result is the same original signal as shown in the red circles below. This all works great so far and proves my wavelength, phases, ampltudes, and equations are all correct. I now want to use the equations for each individual frequency curve to predict a future point. I can do this for every individual equation with no issues as shown below. This is a single frequency from the spectrum. I then sum these predictions, the same as summing the equations, and this is where issues arise. This works only if my wavelength is divisble by my timesteps. So, if I have a range (1, 100) dt=1, and a wavelength of 20, this works. If I change it to (1, 110) it is completely upside down, however everything else still works. What is going on here? It is a LOT of code but I can show anything upon request if needed. Here is the summed signals working.. range (1, 100) wavelength=20 And here it is when I change the timesteps to not be divisible by the wavelengths... range (1, 90) wavelength=20 range (1, 95) wavelength=20 Changing the wavelength would have the same effect as well. What is going on here???

Furthermore, I can show that each individual frequency is predicting correctly. The solid lines left of the prediction, when summed, equal the actual curve. However, the solid lines for the prediction, which are correct as shown, when summed, do not match the actual validation set.

from numpy import pi, sin, arange, mean, abs as absolute, where, nanmin, nanmax, angle, arctan2, sqrt, array, random, append
from matplotlib import pyplot as plt
from scipy import fftpack
import pandas as pd

class GenericSignals():

def __init__(self, time_vec):
self.time_vec = time_vec

def createSingleSignal(self, w=.2, a=10, ph=0):
# x is x axis time vector, w is wavelength, a is amplitude, ph is phase
x = self.time_vec
return a * sin((2 * pi / w * x) + ph)

def combineNSignals(self, *args):
return sum(args)

def getFrequencies(self, signal, dt):
sig_fft = fftpack.fft(signal)
sample_freq = fftpack.fftfreq(len(signal), d=dt)
power = absolute(sig_fft) * dt
return power, sample_freq

def findPeakFrequencies(self, sample_freq, power, n=1):
#n is number of frequencies to get
pos_power_sorted = sorted(pos_power)
dominent = pos_power_sorted[-n:]
# need to reorganize dominent to have proper sorting now
dominant_freq_indices = []
peak_freqs = []
for i in range(len(dominent)):
dominant_freq_indices.append(where(pos_power == dominent[i])[0][0])
peak_freqs.append(freqs[dominant_freq_indices[i]])
return peak_freqs

def createIndividualSignalsForEachFreq(self, signal, sample_freq):
filtered_signals = []
for i in range(len(sample_freq)):
high_freq_fft = fftpack.fft(signal)
high_freq_fft[absolute(sample_freq) < nanmin(sample_freq[i])] = 0
high_freq_fft[absolute(sample_freq) > nanmax(sample_freq[i])] = 0
filtered_sig = fftpack.ifft(high_freq_fft)
filtered_sig -= mean(filtered_sig)
filtered_signals.append(filtered_sig)
return filtered_signals

def getPhaseAmplitudeWavelength(self, signal, freq, sample_freq):
# if statement resolves a divide by 0 runtime warning
if freq == 0:
fixed_freq = .0000000000001
wavelength = 1 / fixed_freq
else:
wavelength = 1 / freq
sig_size = len(signal)
sig_fft = fftpack.fft(signal)
sample_index = where(sample_freq==freq)
phase = (arctan2(sig_fft[sample_index].imag, sig_fft[sample_index].real))[0]
ph = phase + pi/2
amp = (sqrt((sig_fft[sample_index].real * sig_fft[sample_index].real) + (sig_fft[sample_index].imag * sig_fft[sample_index].imag)) / (sig_size / 2))[0]
return ph, amp, wavelength

def getAllPhaseAmplitudeWavelengths(self, all_signals, sample_freq):
wavelengths = []
phases = []
amplitudes = []
i = 0
for individual_signal in all_signals:
phase, amplitude, wavelength = self.getPhaseAmplitudeWavelength(individual_signal, sample_freq[i], sample_freq)
i += 1
wavelengths.append(wavelength)
phases.append(phase)
amplitudes.append(amplitude)
return wavelengths, phases, amplitudes

def eqn(self, signal, wavelength, time_vec, phase, amp):
signals_mean = absolute(mean(signal))
return (amp * sin((2 * pi / wavelength * time_vec) + phase)) + signals_mean

def getEquations(self, wavelength, time_vec, ph, amp):
equations = []
for i in range(len(wavelength)):
equation = (amp[i] * sin((2 * pi / wavelength[i] * time_vec) + ph[i]))
equations.append(equation)
return equations

def predictFuture(self, new_time_vec, equations, wavelength, ph, amp):
# addidtional_step = new_time_vec[-1] + 1
pred = []
for i in range(len(equations)):
pred.append(self.eqn(equations[i], wavelength[i], new_time_vec, ph[i], amp[i]))
return pred

class PreProcessData():

def __init__(self, data):
self.data = data

# separate data into test and validate sets for x and y and for each interval
def createTestValidateSets(self, data, interval_to_predict):
# interval to predict must be less than 1/3 size od dataset
try:
data_x = data[0]
data_y = data[1]
if interval_to_predict * 3 >= len(data_x):
interval_to_predict = 1
if interval_to_predict == 0:
test_x = data_x
val_x = []
test_y = data_y
val_y = []
else:
test_x = data_x[:-interval_to_predict]
val_x = data_x[-interval_to_predict:]
test_y = data_y[:-interval_to_predict]
val_y = data_y[-interval_to_predict:]
return test_x, test_y, val_x, val_y
except Exception as e:
return 'PreProcessData.createTestValidateSets failed: ' + e

if __name__ == "__main__":

dt = 1
time_vec = arange(0, 100, dt)
genSig = GenericSignals(time_vec)
testing_signal = genSig.createSingleSignal(w=20, a=15, ph=0)
testing_signal2 = genSig.createSingleSignal(w=10, a=15, ph=0)
test_signal = testing_signal
predict_interval = 10
preProcessedData = PreProcessData(test_signal)
test_x, test_y, val_x, val_y = preProcessedData.createTestValidateSets([time_vec, test_signal], predict_interval)
power, sample_freq = genSig.getFrequencies(test_y, dt)
individual_signals = genSig.createIndividualSignalsForEachFreq(test_y, sample_freq)
test_index = 5
fft = fftpack.fft(individual_signals[test_index])
fft_phase = arctan2(fft[test_index].imag, fft[test_index].real)

wavelengths, phases, amplitudes = genSig.getAllPhaseAmplitudeWavelengths(individual_signals, sample_freq)
equations = genSig.getEquations(wavelengths, test_x, phases, amplitudes)
index2 = int(len(test_x) / 2 - 1)
all_equations = sum(equations[1:index2])
y_shift = mean(test_y) - mean(all_equations)
all_equations += y_shift
predicted_curves = genSig.predictFuture(val_x, equations, wavelengths, phases, amplitudes)
total_pred = sum(predicted_curves)
total_pred += y_shift
fig, axs = plt.subplots(1, 1, figsize=(15,5))
axs.plot(test_x, individual_signals[test_index], 'r--', label='single exact curve')
axs.plot(test_x, equations[test_index], 'ro', label='single equation curve')
axs.plot(val_x, predicted_curves[test_index], 'r', label='single predcted curve')
# axs.plot(time_vec, test_signal, 'b--', label='test set')
# axs.plot(test_x, all_equations, 'bo', label='ALL equations set')
# axs.plot(val_x, val_y, 'g--', label='val set')
# for i in range(1, 8):
#     axs.plot(val_x, predicted_curves[i], label=f'pred set: {i}')
# axs.plot(val_x, sum(predicted_curves[1:3]), label='pred curves')
axs.plot(val_x, total_pred, 'go', label='predicted set')
axs.legend()
fig.tight_layout()
plt.show()



To clarify, the FFT is working inside the bounds of the FFT. It is outsde that is not working. Below shows the bounds for which the FFT was performed. The stuff to the right is being calculated and summed like a sum of sin waves. This is where the problem exists.

• I'd bet it's just a coding mistake, such as ignoring the phase in the extrapolation. Without seeing the code, it's hard to tell, though.
– MBaz
Apr 6 at 18:09
• I reduced the code down to only the things being used for the most part. Let me know if you have any questions. You can play around with the time_vec and/or test_signal and plot stuff to see what I am talking about. Apr 6 at 18:25
• I added an image that shows my individual frequencies are predicting correctly as well. Very strange. Apr 6 at 18:36
• @MBaz have you had a chance to look at this? This might just be an artifact of how the physics works? At the edges of the FFT things break and therefor you cannot calculate outside of the known data? Apr 7 at 17:14
• I don't have the time right now, sorry -- it's a lot of code to read in a language I don't normally use. Hopefully someone else will chime in. It would help if you could provide the smallest amount of code that reproduces the problem.
– MBaz
Apr 7 at 17:27

I think there's a simple explanation: the expansion of a DFT is periodic. In other words, adding all the sinusoids found by the DFT over a long period of time will produce repetitions of the data given to the DFT. Here's an example that reproduces your situation. The blue dots are the samples given to the DFT, which were taken over a 0.25 second period. The expansion of the sinusoids calculated by the DFT, plotted in red, is also periodic with period 0.25.

Your program works as you expected when an integer number of periods fits exactly into the sampling interval, since no discontinuity is introduced in that case.

Is it possible, in general, to extrapolate a signal using the DFT? I think the best you can do, when a signal is known to consist of a small number $$N$$ of sinusoids, is to keep only the $$N$$ largest DFT bins and extend them -- but the result will not be exact, because a fraction of the input energy will spill to bins that are being ignored. Here's an example with a single sinusoid:

• Thanks for your response! I think we are doing different things though. You are reproducing the signal inside the FFT window, which I have also done successfully. I am trying to continue that signal outside of the FFT range. It is here that I have issues. Apr 13 at 19:27
• Also, unfortunately your last suggestion will not work. When the wavelength is exactly half of the timestep intervals you get 2 dominent waves inverse of one another and the small ones are what matter. for example using 20 wavelength and x=range[0,90]. This produces dominant wavelengths at 22.5 and 18 and everything else gets smaller. The result is 100% out of phase prediction. Apr 13 at 19:32
• Well, in my examples the FFT is taken for samples taken from 0 to 0.25 s -- the red line is an extension outside that range. Anyway, if you ever have a chance to reduce your code or to describe the math behind your algorithm, give me a shout and I'll look again.
– MBaz
Apr 13 at 21:01
• Ok, will do. Can you show your code? I might be able to see more clearly this way also. Apr 13 at 21:15
• It looks like you are showing the exact same problem as me but you have chosen less noticeable terms. Can you try a range from 0 to 90 with dt=1, wavelength=20, and phase=0 for the curve. This will put it completely out of phase for the prediction I think. If you change this wavelength to 30, it will be spot on. Or conversely if you change the range from 0 to 100 it will be spot on as well. Apr 13 at 21:25