Implémentation de l'équation de Wikipedia pour la DFT

10

J'écrivais une implémentation de transformation de Fourier simple et j'ai regardé l' équation DFT sur wikipedia pour référence , quand j'ai remarqué que je faisais quelque chose différemment, et après y avoir pensé, j'ai senti que la version wikipedia devait être incorrecte car il est très simple de penser à un signal que lorsque transformée de Fourier (avec cette équation) retournera un spectre incorrect: parce que l'équation enveloppe le signal autour du plan complexe une seule fois (en raison du avec ), tout signal qui est périodique un nombre pair de fois (tout en enveloppant le plan complexe) n'aura pas de spectre car les pics habituels (tout en faisant le tour du cercle unitaire) qui apparaîtraient pendant un DFT s'annuleraient mutuellement (lorsqu'un nombre pair d'entre eux apparaîtrait). $n/N$ $0<n<N-1$

Pour vérifier cela, j'ai écrit du code qui a produit l'image suivante, qui semble confirmer ce que je pense.

"Le temps en utilisant l'équation" utilise l'équation avec un vecteur de temps (donc le temps auquel été échantillonné par exemple). Il peut être trouvé dans la fonction ci-dessous.

X_{f} = \sum_{n = 0}^{N - 1} x_{n} (\cos (2 π f t_{n}) - i \sin (2 π f t_{n}))

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

t

$t$

t_{n}

$t_n$

x_{n}

$x_n$ ft

L'équation wikipedia, liée ci-dessus, est copiée ici pour référence: Il peut être trouvé dans la fonction .

X_{f} = \sum_{n = 0}^{N - 1} x_{n} (\cos (2 π f \frac{n}{N}) - i \sin (2 π f \frac{n}{N}))

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

De toute évidence, il semble plutôt improbable que j'aurais trouvé au hasard une erreur sur une page wiki aussi prestigieuse. Mais je ne vois pas d'erreur dans ce que j'ai fait?

fourier-transform Nimitz14
la source

Pour mieux comprendre la signification d'un DFT, je vous recommande de lire mes deux premiers articles de blog: "La nature exponentielle du cercle d'unités complexes" ( dsprelated.com/showarticle/754.php ) et "DFT Graphical Interpretation: Centroids des racines pondérées de l'unité "( dsprelated.com/showarticle/768.php ).

Cedron Dawg

Merci je vais jeter un oeil. Je suis honnêtement très surpris de l'attention que cela a eue quand tout cela est dû à un bug très stupide dans mon code.

Nimitz14

Je suis aussi surpris. La chose continue vs discrète est un gros problème cependant. Mon blog est tout au sujet du cas discret sans aucune référence au cas continu qui est différent de l'enseignement du cas discret comme une version échantillonnée du cas continu.

Cedron Dawg

16

Vous avez un bug ft2. Vous incrémentez i, et freqensemble. Ce n'est pas ainsi que vous souhaitez que votre sommation fonctionne. Je me suis amusé à le réparer, mais c'est devenu compliqué. J'ai décidé de le réécrire dans une perspective discrète au lieu d'essayer d'utiliser la terminologie continue. Dans le DFT, le taux d'échantillonnage n'est pas pertinent. Ce qui compte, c'est le nombre d'échantillons utilisés ( N). Les nombres de cases ( k) correspondent alors à la fréquence en unités de cycles par trame. J'ai essayé de vous garder le code aussi intact que possible afin qu'il reste facilement compréhensible pour vous. J'ai également déroulé les boucles de calcul DFT pour, espérons-le, révéler un peu mieux leur nature.

J'espère que cela t'aides.

Ced

import numpy as np
import matplotlib.pyplot as plt 

def ft (t, s, fs):
    freq_step = fs / len (s)
    freqs = np.arange (0, fs / 2, freq_step)
    S = []
    pour freq en freqs:
        réel = np.sum (s * np.cos (2 * np.pi * freq * t)) 
        compl = np.sum (- s * np.sin (2 * np.pi * freq * t)) 
        tmpsum = (réel ** 2 + compl ** 2) ** 0,5 
        S.append (tmpsum)
    return S, freqs

def ft3 (s, N): # Forme d'équation wikipedia plus efficace

    S = []

    tranche = 0,0
    ruban = 2 * np.pi / flotteur (N) 

    pour k dans la plage (N / 2):

        sum_real = 0.0    
        sum_imag = 0.0
        angle = 0,0
        pour n dans la plage (N):
            sum_real + = s [n] * np.cos (angle)
            sum_imag + = -s [n] * np.sin (angle)
            angle + = tranche

        tranche + = ruban
        tmpsum = (sum_real ** 2 + sum_imag ** 2) ** 0,5 
        S.append (tmpsum)

    Retour

def ft4 (s, N): # Utilisation de l'équation wikipedia

    S = []

    pour k dans la plage (N / 2):

        sum_real = 0.0    
        sum_imag = 0.0
        pour n dans la plage (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)

    Retour

def ft5 (s, N): # somme pondérée des racines de l'unité

    ruban = 2 * np.pi / flotteur (N) 

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

    angle = 0,0
    pour r dans la plage (N):
        root_real [r] = np.cos (angle)
        root_imag [r] = -np.sin (angle)
        angle + = ruban

    S = []

    pour k dans la plage (N / 2):

        sum_real = 0.0    
        sum_imag = 0.0
        r = 0

        pour n dans la plage (N):
            sum_real + = s [n] * root_real [r]
            sum_imag + = s [n] * root_imag [r]
            r + = k
            si r> = N: r - = N

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

    Retour

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 ('Signal dans le domaine temporel')
    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 ('Temps utilisant l'équation')
    ax.set_xlabel ('fréquence')
    ax.plot (freqs, S)

    S = pi3 (y, N) 
    ax = fig.add_subplot (313)
    ax.set_title ('Utilisation de l'équation Wiki')
    ax.set_xlabel ('fréquence')
    ax.set_xticks (np.arange (0, freqs [-1], 2)) 
    imprimer len (S), len (freqs)
    ax.plot (freqs, S)

    plt.tight_layout ()
    plt.show ()

principale()

Cedron Dawg
la source

btw vous aviez probablement des problèmes parce que mon code supposait python3;)

Nimitz14

1

@ Nimitz14, Pas un gros problème. J'ai ajouté le "float ()" et un tas de ".0" sur les chiffres. Votre code a bien fonctionné, la seule chose que j'ai dû supprimer était l'instruction "plt.style.use ('ggplot')".

Cedron Dawg

1

@ Nimitz14, j'ai oublié de mentionner, j'ai ajouté une routine ft5 au code qui pré-calcule les racines des valeurs d'unité et montre vraiment comment la DFT est calculée en utilisant les mêmes racines pour chaque bin.

Cedron Dawg

4

je ne vais pas regarder dans votre code. la page de wikipedia semble correct, mais il est un bon exemple de la « guerre des formats » ou « guerre de notation » ou « guerre de style » entre mathématiciens et ingénieurs électriciens. certains, je pense que les maths ont raison. Les EE n'auraient jamais dû adopter " " pour l'unité imaginaire. cela dit, c'est une meilleure expression de la DFT et l'inverse est: $j$

DFT:

X [k] = \sum_{n = 0}^{N - 1} x [n] e^{- j 2 π n k / N}

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

iDFT:

x [n] = \frac{1}{N} \sum_{k = 0}^{N - 1} X [k] e^{j 2 π n k / N}

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

parce que les ingénieurs électriciens qui utilisent DSP aiment utiliser comme séquence d'échantillons en "temps" et comme séquence d'échantillons discrets en "fréquence". les mathématiciens aimeraient mieux ceci: $x[n]$ $X[k]$

DFT:

X_{k} = \sum_{n = 0}^{N - 1} x_{n} e^{- i 2 π n k / N}

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

iDFT:

x_{n} = \frac{1}{N} \sum_{k = 0}^{N - 1} X_{k} e^{i 2 π n k / N}

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

et c'est la même chose que la page wikipedia.

vous devrez peut-être accorder plus d'attention à l'utilisation de ou $+$ $-$ dans l'exposant et comment cela se traduit par $+$ ou $-$ contre la $\sin(\cdot)$ terme.

robert bristow-johnson
la source

3

Si nous utilisions i au lieu de j, nous ne pourrions pas dire ELI l'homme ICE. ELJ l'homme JCE n'a pas la même bague. La civilisation serait menacée

1

Elijah l'homme jus?

robert bristow-johnson

@ user28715 Eh bien, je dans ce cas est en cours pas la racine carrée de moins 1 .... youtube.com/watch?v=2yqjMiFUMlA

Peter K.

0

J'y suis revenu et j'ai essayé de dériver la version discrète qui a aidé à donner plus de sens aux choses:

En quelque sorte $f_k t_n = f(n, k, N)$

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

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

Donc

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

Terminé!

Nimitz14
la source

Implémentation de l'équation de Wikipedia pour la DFT

Réponses: