Implementasi Persamaan Wikipedia untuk DFT

Saya sedang menulis implementasi transformasi fourier sederhana dan melihat persamaan DFT di wikipedia untuk referensi , ketika saya perhatikan bahwa saya melakukan sesuatu yang berbeda, dan setelah memikirkannya merasa bahwa versi wikipedia pasti salah karena sangat mudah untuk memikirkan sebuah sinyal bahwa ketika fourier ditransformasikan (dengan persamaan itu) akan mengembalikan spektrum yang salah: Karena persamaan membungkus sinyal di sekitar bidang kompleks hanya sekali (karena dengan ), setiap sinyal yang periodik beberapa kali genap (saat membungkus bidang kompleks) tidak akan memiliki spektrum seperti puncak biasa (saat mengelilingi lingkaran unit) yang akan muncul selama DFT akan membatalkan satu sama lain (ketika angka genap dari mereka muncul).$n/N$$0<n<N-1$

Untuk memeriksa ini saya menulis beberapa kode yang menghasilkan gambar berikut, yang sepertinya mengkonfirmasi apa yang saya pikirkan.

"Waktu menggunakan persamaan" menggunakan persamaan dengan vektor waktu (jadi waktu di mana disampel misalnya). Ini dapat ditemukan pada fungsi di bawah ini.

$$X_f = \sum_{n=0}^{N-1} x_n (\cos(2\pi f t_n) - i \sin(2\pi f t_n))$$ $t$$t_n$$x_n$ft

Persamaan wikipedia, ditautkan di atas, disalin di sini untuk referensi: Dapat ditemukan dalam fungsi .

$$X_f = \sum_{n=0}^{N-1} x_n \left(\cos\left(2\pi f \frac{n}{N}\right) - i\sin\left(2\pi f \frac{n}{N}\right)\right)$$ ft2

import numpy as np
import matplotlib.pyplot as plt 
plt.style.use('ggplot')

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t)) 
        compl = np.sum(- s * np.sin(2*np.pi*freq * t)) 
        tmpsum = (real**2 + compl**2) ** 0.5 
        S.append(tmpsum)
    return S, freqs

def ft2(s, fs):  # Using wikipedia equation
    nump=len(s)
    freq_step = fs / nump
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for i, freq in enumerate(freqs):
        real = np.sum(s * np.cos(2*np.pi*freq * i/nump))
        compl = np.sum(- s * np.sin(2*np.pi*freq * i/nump))
        tmpsum = (real**2 + compl**2) ** 0.5 
        S.append(tmpsum)
    return S, freqs


def main():
    f = 5 
    fs = 100 
    t = np.linspace(0, 2, 200)
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs) 

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S, freqs = ft2(y, fs) 
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()

Jelas sepertinya agak tidak mungkin bahwa saya akan secara acak menemukan kesalahan pada halaman wiki profil tinggi. Tetapi saya tidak dapat melihat kesalahan dalam apa yang telah saya lakukan?

Anda memiliki bug di ft2. Anda bertambah i, dan freqbersama - sama. Bukan itu yang Anda inginkan agar penjumlahan Anda bekerja. Saya bermain-main dengan memperbaikinya, tetapi itu menjadi berantakan. Saya memutuskan untuk menulis ulang dari perspektif diskrit alih-alih mencoba menggunakan terminologi berkelanjutan. Dalam DFT, laju pengambilan sampel tidak relevan. Yang penting adalah berapa banyak sampel yang digunakan ( N). Nomor bin ( k) kemudian sesuai dengan frekuensi dalam satuan siklus per frame. Saya mencoba untuk menjaga kode Anda semaksimal mungkin sehingga akan tetap mudah dipahami oleh Anda. Saya juga membuka loop perhitungan DFT untuk mudah-mudahan mengungkapkan sifat mereka sedikit lebih baik.

Semoga ini membantu.

Ced

impor numpy sebagai np
impor matplotlib.pyplot sebagai plt 

def ft (t, s, fs):
    freq_step = fs / len (s)
    freqs = np.arange (0, fs / 2, freq_step)
    S = []
    untuk freq dalam freqs:
        real = np.sum (s * np.cos (2 * np.pi * freq * t)) 
        compl = np.sum (- s * np.sin (2 * np.pi * freq * t)) 
        tmpsum = (real ** 2 + compl ** 2) ** 0.5 
        S.append (tmpsum)
    return S, freqs

def ft3 (s, N): # Bentuk persamaan wikipedia yang lebih efisien

    S = []

    slice = 0,0
    sliver = 2 * np.pi / float (N) 

    untuk k dalam kisaran (N / 2):

        sum_real = 0,0    
        sum_imag = 0,0
        sudut = 0,0
        untuk n dalam rentang (N):
            sum_real + = s [n] * np.cos (sudut)
            sum_imag + = -s [n] * np.sin (sudut)
            angle + = slice

        slice + = sliver
        tmpsum = (sum_real ** 2 + sum_imag ** 2) ** 0.5 
        S.append (tmpsum)

    kembali S

def ft4 (s, N): # Menggunakan persamaan wikipedia

    S = []

    untuk k dalam kisaran (N / 2):

        sum_real = 0,0    
        sum_imag = 0,0
        untuk n dalam rentang (N):
            sum_real + = s [n] * np.cos (2 * np.pi * k * n / float (N))
            sum_imag + = -s [n] * np.sin (2 * np.pi * k * n / float (N))

        tmpsum = (sum_real ** 2 + sum_imag ** 2) ** 0.5 
        S.append (tmpsum)

    kembali S

def ft5 (s, N): # Roots of Unity Weighted Sum

    sliver = 2 * np.pi / float (N) 

    root_real = np.zeros (N)
    root_imag = np.zeros (N)

    sudut = 0,0
    untuk r dalam kisaran (N):
        root_real [r] = np.cos (sudut)
        root_imag [r] = -np.sin (sudut)
        angle + = sliver

    S = []

    untuk k dalam kisaran (N / 2):

        sum_real = 0,0    
        sum_imag = 0,0
        r = 0

        untuk n dalam rentang (N):
            sum_real + = s [n] * root_real [r]
            sum_imag + = s [n] * root_imag [r]
            r + = k
            jika r> = N: r - = N

        tmpsum = np.sqrt (sum_real * sum_real + sum_imag * sum_imag)
        S.append (tmpsum)

    kembali S

def main ():

    N = 200
    fs = 100.0

    time_step = 1.0 / fs
    t = np.arange (0, N * time_step, time_step)

    f = 5.0
    y = np.sin (2 * np.pi * f * t) + np.cos (2 * np.pi * f * 2 * t)

    fig = plt.figure ()
    ax = fig.add_subplot (311)
    ax.set_title ('Sinyal dalam domain waktu')
    ax.set_xlabel ('t')
    kapak.plot (t, y)

    S, freqs = ft (t, y, fs) 

    ax = fig.add_subplot (312)
    ax.set_xticks (np.arange (0, freqs [-1], 2)) 
    ax.set_title ('Waktu menggunakan persamaan')
    ax.set_xlabel ('frekuensi')
    kap.plot (freqs, S)

    S = ft3 (y, N) 
    ax = fig.add_subplot (313)
    ax.set_title ('Menggunakan persamaan Wiki')
    ax.set_xlabel ('frekuensi')
    ax.set_xticks (np.arange (0, freqs [-1], 2)) 
    print len ​​(S), len (freqs)
    kap.plot (freqs, S)

    plt.tight_layout ()
    plt.show ()

utama()

saya tidak akan melihat kode Anda. halaman wikipedia terlihat oke, tetapi ini adalah contoh yang bagus dari "perang format" atau "perang notasi" atau "perang gaya" antara matematikawan dan insinyur listrik. beberapa di antaranya, saya pikir orang matematika benar. EE seharusnya tidak pernah diadopsi "$j$"untuk unit imajiner. Yang mengatakan, ini adalah ekspresi DFT yang lebih baik dan kebalikannya adalah:

DFT:

$$X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi nk/N}$$

iDFT:

$$x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] \, e^{j 2 \pi nk/N}$$

karena insinyur listrik yang melakukan DSP suka menggunakan $x[n]$ sebagai urutan sampel dalam "waktu" dan $X[k]$sebagai urutan sampel diskrit dalam "frekuensi". ahli matematika mungkin lebih suka ini:

DFT:

$$X_k = \sum\limits_{n=0}^{N-1} x_n \, e^{-i 2 \pi nk/N}$$

iDFT:

$$x_n = \frac{1}{N} \sum\limits_{k=0}^{N-1} X_k \, e^{i 2 \pi nk/N}$$

dan itu sama dengan halaman wikipedia.

Anda mungkin perlu lebih memperhatikan penggunaan $+$ atau $-$ dalam eksponen dan bagaimana artinya $+$ atau $-$ melawan $\sin(\cdot)$ istilah.

Saya kembali ke ini dan mencoba menurunkan versi diskrit yang membantu membuat segalanya lebih masuk akal:

Entah bagaimana $f_k t_n = f(n, k, N)$

$f_k = \frac{f_s}{N}k$ dan $t_n = \frac{T}{N}n$

$f_s = \frac{N}{T}$

Begitu

$f_k t_n = \frac{f_s}{N}k\frac{T}{N}n = \frac{N}{TN}k\frac{T}{N}n = \frac{kn}{N}$

Selesai!