Аппроксимация волнами sin с использованием DFT на python. Что не так?

Question

Я пишу пропрограмму на python, которая может приближать временные ряды по грехам. Программа использует DFT для поиска волн sin, после чего выбирает грешные волны с наибольшими амплитудами.

Вот мой код:

__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() 

Но у меня есть странный результат:

В программе я создал временную ряд из двух волн греха:

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

И моя программа обнаружила неправильные параметры волн sin:

программа нашла следующие 2 грех волны: [(АМФ = 0,9998029885151699, FRQ = 0,10053096491487339, фаза = 1,1411803525843616), (АМФ = 0.24800925225626422, Frq = 0.40212385965949354, фаза = 0,346757128184013)]

Я suppuse , что моя проблема - неправильное масштабирование волновых параметров, но я не уверен. Есть два места, где программа масштабируется. Первое место создания волн:

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

И второе место sclaling по оси х:

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

Но мой, может быть неправильным. Я пытался изменить масштабирование много раз, и он не решил мою проблему.

Что может быть неправильным с кодом?

Answer 1

Таким образом, эти строки я считаю, не правы:

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])) 

Я считаю, что это должно быть деления на т вместо п, так как вы работаете только с вычислением половину очков. Это устранит проблему амплитуды. Кроме того, при вычислении imag[w] отсутствует отрицательный знак. Принимая во внимание ATAN2 исправить, это будет выглядеть так:

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])) 

Следующий здесь:

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

Дележ по m не прав. amp имеет только половину своих очков, поэтому использование длины усилителя здесь не так.Оно должно быть:

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

Наконец, ваша реконструкция модели есть проблема:

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

Он должен использовать косинус вместо (синусоида вводит сдвиг фазы):

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

Когда я фиксирую все это, модель и ДПФ имеют смысл: