2
$\begingroup$

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. enter image description here 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 enter image description here And here it is when I change the timesteps to not be divisible by the wavelengths... range (1, 90) wavelength=20 enter image description here range (1, 95) wavelength=20 enter image description here 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. enter image description here

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_mask = where(sample_freq > 0)
        freqs = sample_freq[pos_mask]
        pos_power = power[pos_mask]
        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
        # new_time_vec = append(new_time_vec, addidtional_step)[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. enter image description here

$\endgroup$
11
  • $\begingroup$ 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. $\endgroup$
    – MBaz
    Apr 6 at 18:09
  • $\begingroup$ 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. $\endgroup$
    – Justin O
    Apr 6 at 18:25
  • $\begingroup$ I added an image that shows my individual frequencies are predicting correctly as well. Very strange. $\endgroup$
    – Justin O
    Apr 6 at 18:36
  • $\begingroup$ @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? $\endgroup$
    – Justin O
    Apr 7 at 17:14
  • $\begingroup$ 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. $\endgroup$
    – MBaz
    Apr 7 at 17:27

1 Answer 1

1
$\begingroup$

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.

enter image description here

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:

enter image description here

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Justin O
    Apr 13 at 19:27
  • $\begingroup$ 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. $\endgroup$
    – Justin O
    Apr 13 at 19:32
  • $\begingroup$ 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. $\endgroup$
    – MBaz
    Apr 13 at 21:01
  • $\begingroup$ Ok, will do. Can you show your code? I might be able to see more clearly this way also. $\endgroup$
    – Justin O
    Apr 13 at 21:15
  • $\begingroup$ 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. $\endgroup$
    – Justin O
    Apr 13 at 21:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.