Giocando con la DFT

[IlGuru]

Ci sono alcune cose che riempiono di meraviglia la mia mente ignorante ogni volta che torno a giocarci, una è l'insieme di Mandelbrot, l'altra è la trasformata di Fourier.

Per potermici trastullare, ho scritto una classe python che calcola la TF discreta e la sua inversa prelevando i dati da una lista di numeri complessi, ed una classe per disegnare dei semplici grafici, che non disegna ancora i numeri quindi per quanto riguarda tutte le statistiche vi dovete fidare dei dati mostrati nei log :)

Per ciascuna funzione ho disegnato tre grafici, uno in disegno il campionamento della funzione trasformata ( in viola ) , uno in cui disegno componenti reali e complesse della dft ( in giallo ), e l'ultimo in cui disegno le componenti raggio e fase della dft in forma polare ( in azzurro ).

Con mia meraviglia, vedo che nel dominio delle frequenze oltre alle componenti che mi aspetto, in fondo al grafico delle frequenze ho due picchi di cui non mi spiego la provenienza. E' palese che abbiamo a che fare con il numero di campioni, ovvero alla frequenza di campionamento del mio frame.

Es:
Numero di campioni: 50




Moltiplicazione di due funzioni seno, una a frequenza 2 ed una a frequenza 10
Qui correttamente trovo un picco a frequenza 8 ed uno a frequenza 12, ma anche uno a frequenza 38 ed uno a frequenza 42:

sen_mul.png (4.04 KiB) Osservato 2646 volte

Qui la parte interessante del log:

Codice: Seleziona tutto

- Calcolo coordinate polari array dft
F:    8 - Rag:  12.5000 Phi:  -0.0000
F:   12 - Rag:  12.5000 Phi:   3.1416
F:   38 - Rag:  12.5000 Phi:   3.1416
F:   42 - Rag:  12.5000 Phi:  -0.0000


Altro esempio:

Come prima, ma somma di due funzioni seno. Qui trovo un picco a frequenza 2, uno a frequenza 10, uno a 40 ed uno a 48

sen_sum.png (3.86 KiB) Osservato 2646 volte

Codice: Seleziona tutto

- Calcolo coordinate polari array dft
F:    2 - Rag:  12.5000 Phi:  -1.5708
F:   10 - Rag:  12.5000 Phi:  -1.5708
F:   40 - Rag:  12.5000 Phi:   1.5708
F:   48 - Rag:  12.5000 Phi:   1.5708


Come si spiegano quei due picchi? E' un bug della mia classe?

Il codice, per quanto buggato, si trova qui:

Codice: Seleziona tutto

git clone https://StefanoBusnelli@bitbucket.org/StefanoBusnelli/pychart.git
git clone https://StefanoBusnelli@bitbucket.org/StefanoBusnelli/pyfourier.git


[xyz]

Stai calcolando una DFT di un segnale reale quindi il segnale trasformato deve essere pari, questo si traduce nella DFT in una serie di valori simmetrici (per la precisione simmetria Hermitiana) presenti all'inizio del vettore (array[i]) e alla fine (array[N - i]). I dettagli matematici sono spiegati qui:

https://dsp.stackexchange.com/questions ... nswer-4827

Il secondo repository GIT indicato prima è privato.

Io per eseguire questo tipo di calcolo in Python utilizzo le seguenti librerie molto diffuse e collaudate:

https://numpy.org
https://www.scipy.org
https://matplotlib.org

Numpy ha due funzioni importanti per mettere in ordine i dati della FFT (algoritmo veloce per il calcolo di una DFT): fftfreq crea l'array con le frequenze e fftshift riordina l'array della FFT spostando i dati in modo da renderli simmetrici al centro.

Questo è un esempio di codice (Scipy in questo caso non è utilizzato):

Codice: Seleziona tutto

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt


def main():

    samples = 50
    t_start = 0
    t_end = 1

    # time domain

    t = np.linspace(t_start, t_end, samples)
    y = np.sin(2 * (2 * np.pi * t)) * np.sin(10 * (2 * np.pi * t))

    # frequency domain

    fft_x = np.fft.fftshift(np.fft.fftfreq(len(y)))
    fft_y = np.fft.fftshift(np.fft.fft(y))

    (fft_min, fft_max) = (fft_x[0], fft_x[-1])

    # plot

    fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1)
    fig.subplots_adjust(hspace=1.0)

    ax1.plot(t, y, 'go-')
    ax1.set_xlim(t_start, t_end)
    ax1.set_title('Time domain')
    ax1.set_xlabel('time')
    ax1.set_ylabel('f(t)')
    ax1.grid(True)

    ax2.plot(fft_x, np.real(fft_y), 'b', label='Re')
    ax2.plot(fft_x, np.imag(fft_y), 'r:', label='Im')
    ax2.plot(fft_x, np.abs(fft_y), 'k--', label='abs')
    ax2.set_xlim(fft_min, fft_max)
    ax2.set_title('Frequency domain')
    ax2.set_xlabel('frequency')
    ax2.set_ylabel('F(f)')
    ax2.grid(True)
    ax2.legend()

    ax3.stem(fft_x, np.real(fft_y))
    ax3.set_xlim(fft_min, fft_max)
    ax3.set_title('Real')
    ax3.set_xlabel('frequency')
    ax3.set_ylabel('Re(F(f))')
    ax3.grid(True)

    ax4.stem(fft_x, np.imag(fft_y))
    ax4.set_xlim(fft_min, fft_max)
    ax4.set_title('Imaginary')
    ax4.set_xlabel('frequency')
    ax4.set_ylabel('Im(F(f))')
    ax4.grid(True)

    ax5.stem(fft_x, np.imag(fft_y))
    ax5.set_xlim(fft_min, fft_max)
    ax5.set_title('Absolute')
    ax5.set_xlabel('frequency')
    ax5.set_ylabel('abs(F(f))')
    ax5.grid(True)

    #plt.savefig("plot.png")
    plt.show()


if __name__ == '__main__':
    main()

Plot:

FFT

[IlGuru]

Grazie per il link xyz, mi hai fornito ottimo materiale su cui approfondire.
Non mi ero accorto che il repository fosse privato, ho rimosso il flag.