Approximation by sin waves using DFT on python. What's wrong?

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:

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?

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: