Approximation par ondes sinus à l'aide de DFT sur python. Qu'est-ce qui ne va pas?

Question

J'écris le prorgramme sur python qui peut approximer des séries temporelles par des ondes sinueuses. Le programme utilise DFT pour trouver les ondes sinusoïdales, après quoi il choisit les ondes sinusoïdales avec les plus grandes amplitudes.

Voici mon 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() 

Mais j'ai un résultat étrange:

Dans le programme que j'ai créé une série chronologique de deux ondes péché:

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

Et mon programme a trouvé les mauvais paramètres des ondes sinueuses:

le programme trouve les 2 prochaines ondes sin: [(amp = 0,9998029885151699, FRQ = ,10053096491487339 phase = 1,1411803525843616), (ampères = 0,24800925225626422, FRQ = 0,40212385965949354 phase = ,346757128184013)]

I suppuse , que mon problème est une mauvaise mise à l'échelle des paramètres d'onde, mais je ne suis pas sûr. Il y a deux endroits où le programme effectue la mise à l'échelle. La première place crée des vagues:

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

Et la deuxième place est sclaling de l'axe des x:

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

Mais mon suppose peut-être tort. J'ai essayé de changer de taille plusieurs fois, et cela n'a pas résolu mon problème.

Qu'est-ce qui ne va pas avec le code?

Answer 1

Ainsi, ces lignes, je crois, sont mal:

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

Je crois qu'il devrait être en divisant par m au lieu de n, puisque vous ne travaillez avec calcul de la moitié des points. Cela va régler le problème d'amplitude. En outre, le calcul de imag[w] manque un signe négatif. Compte tenu de la solution atan2, il ressemblerait à ceci:

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

Le prochain est ici:

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

La division par m est pas juste. amp a seulement la moitié des points qu'il devrait, donc en utilisant la longueur de l'ampli n'est pas juste ici.Il devrait être:

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

Enfin, votre reconstruction de modèle a un problème:

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

Il devrait utiliser à la place cosinus (sinus introduit un décalage de phase):

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

Quand je fixe Tout cela, le modèle et le DFT ont tous deux un sens: