Transformada Discreta de Fourier¶
Implementação da DFT Matricial
Aluno: Alysson Machado de Oliveira Barbosa
In [1]:
# importação dos pacotes utilitários
import matplotlib.pyplot as plt
import numpy as np
import warnings
import time
# desabilita todos os warnings
warnings.filterwarnings("ignore")
Lembrando que se $x(n)$ é um sinal no tempo de discreto e aperiódico, sua Transformada de Fourier é dada por $$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n} $$ e a inversa $$ x(n) = \frac{1}{2\pi}\int_{2\pi}X(e^{j\omega})e^{j\omega n}d\omega $$ $X(e^{j\omega})$ é uma função contínua de $\omega$.
Exemplo:
Seja o sinal $$ x(n) = \left\{ \begin{array}{c c} 1, & 0 \leq n < L \\ 0, & c.c. \end{array} \right. $$ que é limitado no tempo com apenas $L$ valores constantes e não nulos. Calculando a Transformada de Fourier para $L=8$ obtemos $$ X(e^{j\omega}) = \frac{e^{-j8\omega} - 1}{e^{-j\omega} - 1} $$ cujo módulo está representado na Figura abaixo:
In [2]:
def x(w, L):
'''
Representação do sinal x[n] no domínio da frequência
Inputs:
w (int) -> frequência em radianos.
L (int) -> valor máximo de n cujo x[n] = 1, considerando que n > 0.
'''
#
numerator = np.exp(-1j * L * w) - 1
denominator = np.exp(-1j * w) - 1
return numerator / denominator
# definindo os valores no domínio da frequência (radianos)
w = np.linspace(0, 2 * np.pi, 200)
# configurações para a plotagem do gráfico
plt.plot(w, np.abs(x(w, L = 8)))
plt.xlabel('Frequência (radianos)')
plt.ylabel('$|x(e^{jw})|$')
plt.title('Módulo da função no domínio da frequência')
plt.grid(True)
plt.show()
O módulo representado na figura acima é uma função contínua de $\omega$, e para realizar a análise espectral usando um processador digital, devemos amostrá-la, o que nos leva a questão de quantas amostras devem ser obtidas.
Por exemplo, se escolhermos $N = 8$, vamos amostras tomar nas frequências $\omega_k = 2\pi k/N$ e obteremos os valores calculados abaixo.
Podemos escrever $$ X[k] = X(e^{j\omega})|_{\omega=\frac{2\pi k}{N}},\ k=0,1,2,..., N-1 $$ e como $X(e^{j\omega})$ é periódica, basta tomar amostras no intervalo de 0 a $2\pi$.
Estamos considerando que o sinal $x(n)$ tem $L$ valores ($x[n] = 0$ para $n < 0 $ e $N\geq L$), ao tomar $N$ amostas do espectro temos $$ X[k] = X(e^{j\frac{2\pi k}{N}}) = \sum\limits_{n=0}^{N-1}x[n]e^{-j\frac{2\pi k}{N}n} $$ ou ainda como $$ X[k] =\left\{ \begin{array}{c c} \sum\limits_{n=0}^{N-1}x[n]W_N^{kn} & k = 0,1, ..., N-1 \\ 0, & c.c \end{array} \right. $$ e $$ x[n] = \left\{ \begin{array}{c c} \frac{1}{N}\sum\limits_{k=0}^{N-1}X[k]W_{N}^{-kn}, & n = 0,1,...,N-1 \\ 0, & c.c. \end{array}\right. $$ sendo $W_N = e^{-j\frac{2\pi}{N}}$.
Definindo a Transformada Discreta de Fourier¶
In [3]:
def dft_sum(time_function, N):
'''
Define a DFT usando procedimentos de repetição com somatórios
Inputs:
time_function (list) -> lista com os valores discretos x[n] da função.
N (int) -> define o espaçamento dos valores discretos x[k] obtidos com a DFT.
Outputs:
x_freq (list) -> valores da DFT x[k] para x[n].
'''
# inicializa uma lista para adicionar os valores de x[k]
x_freq = np.zeros((N))
# aplica o zero padding em x[n] para os casos em que N > len(x[n])
if len(time_function) < N:
time_function.extend(np.zeros(N - len(time_function)))
# somatório para x[0], ..., x[N-1], 0 <= k < N - 1
for k in range(0, N):
# intera as componentes de cada somatório de x[k], 0 <= n < N - 1
x_k_values = list()
for n in range(0, N):
x_k_values.append(np.multiply(time_function[n], np.exp(-1j * ((2 * np.pi) / N) * (k * n))))
x_freq[k] = np.sum(x_k_values)
return x_freq
In [4]:
def dft_matrix(time_function, N):
'''
Define a DFT usando procedimentos de multiplicação matricial.
Inputs:
time_function (list) -> lista com os valores discretos x[n] da função.
N (int) -> define o espaçamento dos valores discretos x[k] obtidos com a DFT.
Outputs:
x_freq (list) -> valores da DFT x[k] para x[n].
'''
# aplica o zero padding em x[n] para os casos em que N > len(x[n])
if len(time_function) < N:
time_function.extend(np.zeros(N - len(time_function)))
# define todos os valores possíveis de n, correspondentes a quantidade máxima de linhas da matriz
n = np.arange(N)
# define todos os valores possíveis de k, correspondentes a quantidade máxima de colunas da matriz
k = n.reshape((N, 1))
# define a matrix que mapeia os elementos de x[n] para x[k]
W = np.exp(-1j * ((2 * np.pi) / N) * (k * n))
# realiza o mapeamento da função x[n] para o domínio x[k]
x_freq = np.dot(W, time_function)
return x_freq
In [5]:
def visualize_dfts(N):
'''
Função para visualizar as DFT pelo método de somatórios e matricial.
Inputs:
N (int) -> define o espaçamento dos valores discretos x[k] obtidos com a DFT.
'''
# aplica a DFT por somatórios
discrete_signal1 = dft_sum([1.,1.,1.,1.,1.,1.,1.,1.], N)
# aplica a DFT por matrizes
discrete_signal2 = dft_matrix([1.,1.,1.,1.,1.,1.,1.,1.], N)
# criando os objetos para plotar duas imagens
fig, axs = plt.subplots(1, 2 , figsize=(16, 6))
# plotando o sinal x(e^(jw))
axs[0].plot(w, x(w, L = 8), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal1)))
# plota os pontos obtidos com a DFT para o método de somatórios
axs[0].scatter(dft_x_axis, discrete_signal1, color = 'red', marker = 'o', label = 'DFT Somatórios')
# define as legendas, títulos e grades da imagem
axs[0].set_xlabel('Frequência (radianos)')
axs[0].set_ylabel('$X(e^{jw})$')
axs[0].set_title(f'Função no domínio da frequência')
axs[0].grid(True)
axs[0].legend()
# plotando o sinal x(e^(jw))
axs[1].plot(w, x(w, L = 8), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal2)))
# plota os pontos obtidos com a DFT para o método de multiplicação matricial
axs[1].scatter(dft_x_axis, discrete_signal2, color = 'red', marker = 'o', label = 'DFT Matricial')
# define as legendas, títulos e grades da imagem
axs[1].set_xlabel('Frequência (radianos)')
axs[1].set_ylabel('$X(e^{jw})$')
axs[1].set_title(f'Função no domínio da frequência')
axs[1].grid(True)
axs[1].legend()
# adiciona um título geral para a figura
fig.suptitle(f'N = {N}')
# Ajustando o layout para evitar sobreposição
plt.tight_layout()
# Exibindo os gráficos
plt.show()
In [6]:
visualize_dfts(N = 8)
In [7]:
visualize_dfts(N = 32)
In [8]:
visualize_dfts(N = 128)
In [9]:
visualize_dfts(N = 512)
Análise das DFT para um impulso e um impulso deslocado¶
In [10]:
def x_impulse(w):
'''
Representação do sinal x[n] no domínio da frequência.
Inputs:
w (int) -> frequência em radianos.
'''
return np.exp(-1j * w * 0)
def x_impulse_shift(w, no):
'''
Representação do sinal x[n] no domínio da frequência.
Inputs:
w (int) -> frequência em radianos.
no (int) -> deslocamento do impulson.
'''
return np.exp(-1j * w * no)
In [11]:
# criando os objetos para plotar duas imagens
fig, axs = plt.subplots(1, 2 , figsize=(16, 6))
discrete_signal3 = dft_sum([1.], N = 16)
discrete_signal4 = dft_matrix([1.], N = 16)
# definindo os valores no domínio da frequência (radianos)
w = np.linspace(0, 2 * np.pi, 200)
# plotando o módulo do sinal x(e^(jw))
axs[0].plot(w, np.abs(x_impulse(w)), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal3)))
# plota os pontos obtidos com a DFT para o método de somatórios
axs[0].scatter(dft_x_axis, np.abs(discrete_signal4), color = 'red', marker = 'o', label = 'DFT Matricial')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
axs[0].set_xlabel('Frequência (radianos)')
axs[0].set_ylabel('$|X(e^{jw})|$')
axs[0].set_title(f'Módulo da função FFT($\delta(n)$) no domínio da frequência')
axs[0].set_ylim([-2, 2])
axs[0].grid(True)
axs[0].legend()
# plotando o ângulo do sinal x(e^(jw))
axs[1].plot(w, np.angle(x_impulse(w)), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal3)))
# plota os pontos obtidos com a DFT para o método de somatórios
axs[1].scatter(dft_x_axis, np.angle(discrete_signal4), color = 'red', marker = 'o', label = 'DFT Matricial')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
axs[1].set_xlabel('Frequência (radianos)')
axs[1].set_ylabel('$\phi$')
axs[1].set_title(f'Ângulo da função FFT($\delta(n)$) no domínio da frequência')
axs[1].set_ylim([-2, 2])
axs[1].grid(True)
axs[1].legend()
# adiciona um título geral para a figura
fig.suptitle(f'N = {64}')
# Ajustando o layout para evitar sobreposição
plt.tight_layout()
# Exibindo os gráficos
plt.show()
In [12]:
# criando os objetos para plotar duas imagens
fig, axs = plt.subplots(1, 2 , figsize=(16, 6))
discrete_signal3 = dft_sum([0, 1., 0, 0], N = 32)
discrete_signal4 = dft_matrix([0, 1., 0, 0], N = 32)
# definindo os valores no domínio da frequência (radianos)
w = np.linspace(0, 2 * np.pi, 200)
# plotando o módulo do sinal x(e^(jw))
axs[0].plot(w, np.abs(x_impulse_shift(w, no = 1)), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal3)))
# plota os pontos obtidos com a DFT para o método de somatórios
axs[0].scatter(dft_x_axis, np.abs(discrete_signal4), color = 'red', marker = 'o', label = 'DFT Matricial')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
axs[0].set_xlabel('Frequência (radianos)')
axs[0].set_ylabel('$|X(e^{jw})|$')
axs[0].set_title(f'Módulo da função FFT($\delta(n-1)$) no domínio da frequência')
axs[0].grid(True)
axs[0].legend()
# plotando o ângulo do sinal x(e^(jw))
axs[1].plot(w, np.angle(x_impulse_shift(w, no = 1)), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal3)))
# plota os pontos obtidos com a DFT para o método de somatórios
axs[1].scatter(dft_x_axis, np.angle(discrete_signal4), color = 'red', marker = 'o', label = 'DFT Matricial')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
axs[1].set_xlabel('Frequência (radianos)')
axs[1].set_ylabel('$\phi$')
axs[1].set_title(f'Ângulo da função FFT($\delta(n-1)$) no domínio da frequência')
axs[1].grid(True)
axs[1].legend()
# adiciona um título geral para a figura
fig.suptitle(f'N = {32}')
# Ajustando o layout para evitar sobreposição
plt.tight_layout()
# Exibindo os gráficos
plt.show()
In [13]:
# criando os objetos para plotar duas imagens
fig, axs = plt.subplots(1, 2 , figsize=(16, 6))
discrete_signal3 = dft_sum([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1., 0, 0, 0,], N = 128)
discrete_signal4 = dft_matrix([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1., 0, 0, 0], N = 128)
# definindo os valores no domínio da frequência (radianos)
w = np.linspace(0, 2 * np.pi, 200)
# plotando o módulo do sinal x(e^(jw))
axs[0].plot(w, np.abs(x_impulse_shift(w, no = 10)), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal3)))
# plota os pontos obtidos com a DFT para o método de somatórios
axs[0].scatter(dft_x_axis, np.abs(discrete_signal4), color = 'red', marker = 'o', label = 'DFT Matricial')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
axs[0].set_xlabel('Frequência (radianos)')
axs[0].set_ylabel('$|X(e^{jw})|$')
axs[0].set_title(f'Módulo da função FFT($\delta(n-10)$) no domínio da frequência')
axs[0].grid(True)
axs[0].legend()
# plotando o ângulo do sinal x(e^(jw))
axs[1].plot(w, np.angle(x_impulse_shift(w, no = 10)), label = 'FFT')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
dft_x_axis = np.arange(0, np.max(w), (np.max(w) / len(discrete_signal3)))
# plota os pontos obtidos com a DFT para o método de somatórios
axs[1].scatter(dft_x_axis, np.angle(discrete_signal4), color = 'red', marker = 'o', label = 'DFT Matricial')
# define os valores no eixo x, igualmente espaçados para a resolução do sinal x(e^(jw))
axs[1].set_xlabel('Frequência (radianos)')
axs[1].set_ylabel('$\phi$')
axs[1].set_title(f'Ângulo da função FFT($\delta(n-10)$) no domínio da frequência')
axs[1].grid(True)
axs[1].legend()
# adiciona um título geral para a figura
fig.suptitle(f'N = {128}')
# Ajustando o layout para evitar sobreposição
plt.tight_layout()
# Exibindo os gráficos
plt.show()
Gráfico do Tempo de Cálculo para Diferentes Valores de N¶
In [23]:
signal = [1,1,1,1,1,1,1,1]
time_inference = list()
for increase in range(len(signal), 1000):
init_time = time.time()
result = dft_matrix(signal, increase)
final_time = time.time()
time_inference.append((final_time - init_time) * 1000)
In [40]:
x = np.arange(8, len(time_inference) + 8)
y = time_inference
grau = 3
coeficientes = np.polyfit(x, y, grau)
polinomio = np.poly1d(coeficientes)
y_plot = polinomio(x)
# Plotagem dos dados originais e da aproximação polinomial
plt.plot(x, y, 'o', label='Tempo de Inferência para N')
plt.plot(x, y_plot, label=f'Aproximação Polinomial (grau {grau})')
plt.xlabel('N')
plt.ylabel('Tempo (segundos)')
plt.title('Gráfico o tempo de execução em função do tamanho da DFT (N)')
plt.legend()
plt.grid(True)
plt.show()
