在 python 上使用 DFT 通过 sin 波近似。怎么了？

Question

我正在 python 上编写程序，可以通过 sin 波来近似时间序列。该程序使用DFT寻找正弦波，然后选择振幅最大的正弦波。

这是我的代码：

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

而且我的程序发现了错误的正弦波参数：

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

我想，我的问题是波参数的缩放比例错误，但我不确定。程序在两个地方进行缩放。第一名是造波：

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

第二名是x轴的缩放：

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

我认为它应该除以 m 而不是 n，因为您只计算了一半的点数。这将解决振幅问题。此外，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 只有应有的一半点数，因此这里不适合使用 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)

当我解决所有这些问题时，模型和 DFT 都有意义：

在 python 上使用 DFT 通过 sin 波近似。怎么了？

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

python

approximation

dft

trigonometry

continuous-fourier