Approximation by sin waves using DFT on python. What's wrong?
Question: I'm writing the prorgram on python that can approximate time series by sin
waves. The program uses DFT to find sin waves, after that it chooses sin waves
with biggest amplitudes.

Here's my code:

    
    
    __author__ = 'FATVVS'
    
    import math
    
    
    # Wave - (amplitude,frequency,phase)
    # This class was created to sort sin waves:
    # - by anplitude( set freq_sort=False)
    # - by frequency (set freq_sort=True)
    class Wave:
        #flag for choosing sort mode:
        # False-sort by amplitude
        # True-by frequency
        freq_sort = False
    
        def __init__(self, amp, freq, phase):
            self.freq = freq #frequency
            self.amp = amp #amplitude
            self.phase = phase
    
        def __lt__(self, other):
            if self.freq_sort:
                return self.freq < other.freq
            else:
                return self.amp < other.amp
    
        def __gt__(self, other):
            if self.freq_sort:
                return self.freq > other.freq
            else:
                return self.amp > other.amp
    
        def __le__(self, other):
            if self.freq_sort:
                return self.freq <= other.freq
            else:
                return self.amp <= other.amp
    
        def __ge__(self, other):
            if self.freq_sort:
                return self.freq >= other.freq
            else:
                return self.amp >= other.amp
    
        def __str__(self):
            s = "(amp=" + str(self.amp) + ",frq=" + str(self.freq) + ",phase=" + str(self.phase) + ")"
            return s
    
        def __repr__(self):
            return self.__str__()
    
    
    #Discrete Fourier Transform
    def dft(series: list):
        n = len(series)
        m = int(n / 2)
        real = [0 for _ in range(n)]
        imag = [0 for _ in range(n)]
        amplitude = []
        phase = []
        angle_const = 2 * math.pi / n
        for w in range(m):
            a = w * angle_const
            for t in range(n):
                real[w] += series[t] * math.cos(a * t)
                imag[w] += series[t] * math.sin(a * t)
            amplitude.append(math.sqrt(real[w] * real[w] + imag[w] * imag[w]) / n)
            phase.append(math.atan(imag[w] / real[w]))
        return amplitude, phase
    
    
    #extract waves from time series
    # series - time series
    # num - number of waves
    def get_waves(series: list, num):
        amp, phase = dft(series)
        m = len(amp)
        waves = []
        for i in range(m):
            waves.append(Wave(amp[i], 2 * math.pi * i / m, phase[i]))
        waves.sort()
        waves.reverse()
        waves = waves[0:num]#extract best waves
        print("the program found the next %s sin waves:"%(num))
        print(waves)#print best waves
        return waves
    
    #approximation by sin waves
    #series - time series
    #num- number of sin waves
    def sin_waves_appr(series: list, num):
        n = len(series)
        freq = get_waves(series, num)
        m = len(freq)
        model = []
        for i in range(n):
            summ = 0
            for j in range(m): #sum by sin waves
                summ += freq[j].amp * math.sin(freq[j].freq * i + freq[j].phase)
            model.append(summ)
        return model
    
    
    if __name__ == '__main__':
        import matplotlib.pyplot as plt
    
        N = 500  # length of time series
        num = 2  # number of sin wawes, that we want to find
    
        #y - generate time series
        y = [2 * math.sin(0.05 * t + 0.5) + 0.5 * math.sin(0.2 * t + 1.5) for t in range(N)]
        model = sin_waves_appr(y, num) #generate approximation model
    
        ## ------------------plotting-----------------
    
        plt.figure(1)
    
        # plotting of time series and his approximation model
        plt.subplot(211)
        h_signal, = plt.plot(y, label='source timeseries')
        h_model, = plt.plot(model, label='model', linestyle='--')
        plt.legend(handles=[h_signal, h_model])
        plt.grid()
    
        # plotting of spectre
        amp, _ = dft(y)
        xaxis = [2*math.pi*i / N for i in range(len(amp))]
        plt.subplot(212)
        h_freq, = plt.plot(xaxis, amp, label='spectre')
        plt.legend(handles=[h_freq])
        plt.grid()
    
        plt.show()
    

But I've got a strange result: [![enter image description
here](http://i.stack.imgur.com/mIonA.png)](http://i.stack.imgur.com/mIonA.png)

In the program I've created a time series from two sin waves:

    
    
    y = [2 * math.sin(0.05 * t + 0.5) + 0.5 * math.sin(0.2 * t + 1.5) for t in range(N)]
    

And my program found wrong parameters of the sin waves:

> the program found the next 2 sin waves:
> [(amp=0.9998029885151699,frq=0.10053096491487339,phase=1.1411803525843616),
> (amp=0.24800925225626422,frq=0.40212385965949354,phase=0.346757128184013)]

I suppuse, that my problem is wrong scaling of wave parameters, but I'm not
sure. There're two places, where the program does scaling. The first place is
creating of waves:

    
    
    for i in range(m):
        waves.append(Wave(amp[i], 2 * math.pi * i / m, phase[i]))
    

And the second place is sclaling of the x-axis:

    
    
    xaxis = [2*math.pi*i / N for i in range(len(amp))]
    

But my suppose may be wrong. I've tried to change scaling many times, and it
haven't solved my problem.

What may be wrong with the code?


Answer: So, these lines I believe are wrong:

    
    
    for t in range(n):
        real[w] += series[t] * math.cos(a * t)
        imag[w] += series[t] * math.sin(a * t)
    amplitude.append(math.sqrt(real[w] * real[w] + imag[w] * imag[w]) / n)
    phase.append(math.atan(imag[w] / real[w]))
    

I believe it should be dividing by m instead of n, since you are only working
with computing half the points. That will fix the amplitude problem. Also, the
computation of `imag[w]` is missing a negative sign. Taking into account the
atan2 fix, it would look like:

    
    
    for t in range(n):
        real[w] += series[t] * math.cos(a * t)
        imag[w] += -1 * series[t] * math.sin(a * t)
    amplitude.append(math.sqrt(real[w] * real[w] + imag[w] * imag[w]) / m)
    phase.append(math.atan2(imag[w], real[w]))
    

The next one is here:

    
    
    for i in range(m):
        waves.append(Wave(amp[i], 2 * math.pi * i / m, phase[i]))
    

The divide by `m` is not right. `amp` has only half the points it should, so
using the length of amp isn't right here. It should be:

    
    
    for i in range(m):
        waves.append(Wave(amp[i], 2 * math.pi * i / (m * 2), phase[i]))
    

Finally, your model reconstruction has a problem:

    
    
        for j in range(m): #sum by sin waves
            summ += freq[j].amp * math.sin(freq[j].freq * i + freq[j].phase)
    

It should use cosine instead (sine introduces a phase shift):

    
    
        for j in range(m): #sum by cos waves
            summ += freq[j].amp * math.cos(freq[j].freq * i + freq[j].phase)
    

When I fix all of that, the model and the DFT both make sense:

[![Picture of spectrum and time
models](http://i.stack.imgur.com/xXLsG.png)](http://i.stack.imgur.com/xXLsG.png)

