Создание функции для обратного преобразования Фурье

Я пытаюсь выполнить обратное преобразование Фурье, создав свою собственную функцию. Это функция преобразования Фурье моего временного ряда, которая работает нормально.

def DFT(x, frequencies):
    N1 = x.size
    k = frequencies.size
    t = np.arange(N1)
    fourier = np.zeros(k)
    for i in range(0,k):
        fourier[i] = np.dot(x,(1/k)*np.exp(-2j*np.pi*frequencies[i]*t))
    return fourier

Это мой исходный сигнал (просто синусоида):

N2 = 1*10**6
time = np.arange(0,N2,1000)
lam = .1*10**6
f = 1/lam
freq = np.arange(0,.05,.0001)


signal = np.sin(2*np.pi*f*time)

И спектр мощности построен с использованием моего ДПФ (функция Фурье):

plt.plot(freq, np.abs(DFT(signal,freq))**2)
plt.xlabel('Frequency')
plt.title('Spectrum of Sine Wave')
plt.grid()
plt.show()

но когда я пытаюсь применить свою функцию для обратного преобразования Фурье, я не получаю исходную синусоиду:

def IFT(fft, frequencies):
    N = fft.size
    k = frequencies.size
    n = np.arange(N)
    inverse_fourier = np.zeros(k)
    for i in range(0,k):
        inverse_fourier[i] = np.dot(fft,np.exp((-2j*np.pi*frequencies[i]*n)/N)) #[None,:]
    return inverse_fourier

Что не так с моей функцией? Я не получаю ошибок, но возвращаемый сигнал совершенно неверен.

python numpy jupyter fft

19.12.2020 11:58

Почему в Python есть оператор "pass"?

Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.

Некоторые методы, о которых вы не знали, что они существуют в Python

Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...

Основы Python Часть I

Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?

LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа

Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:

Оптимизация кода с помощью тернарного оператора Python

И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это

Советы по эффективной веб-разработке с помощью Python

Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.

641

Перейти к ответу Данный вопрос помечен как решенный

Ответы 2

Ответ принят как подходящий

В вашем обратном показателе степень должна быть положительной (2j вместо -2j) и отбросить /N, что дает (добавлены графики для демонстрации):

import numpy as np
import matplotlib.pyplot as plt

def DFT(x, frequencies):
    N1 = x.size
    k = frequencies.size
    t = np.arange(N1)
    fourier = np.zeros(k)
    for i in range(0,k):
        fourier[i] = np.dot(x, (1/k)*np.exp(-2j*np.pi*frequencies[i]*t))
    return fourier

def IFT(fft, frequencies):
    N = fft.size
    k = frequencies.size
    n = np.arange(N)
    inverse_fourier = np.zeros(k)
    for i in range(0,k):
        inverse_fourier[i] = np.dot(fft, np.exp((2j*np.pi*frequencies[i]*n))) #[None,:]
    return inverse_fourier

N2 = 1*10**6
time = np.arange(0,N2,2000)
lam = .1*10**6
f = 1/lam
freq = np.arange(0,.05,.0001)

signal = np.sin(2*np.pi*f*time)

plt.plot(time, signal)
plt.xlabel('Time')
plt.title('Sine Wave')
plt.grid()
plt.show()

dft = DFT(signal, freq)

plt.plot(freq, np.abs(dft)**2)
plt.xlabel('Frequency')
plt.title('Spectrum of Sine Wave')
plt.grid()
plt.show()

plt.plot(time, IFT(dft, freq))
plt.xlabel('Time')
plt.title('Sine Wave')
plt.grid()
plt.show()

что дает (первый график греха опущен):

19.12.2020 15:11

Запустив код, вы должны получить следующее предупреждение:

ComplexWarning: Casting complex values to real discards the imaginary part
  fourier[i] = np.dot(x,(1/k)*np.exp(-2j*np.pi*frequencies[i]*t))

Поскольку результирующее преобразование Фурье должно быть комплексным, это предупреждение должно вызывать опасения. Чтобы избавиться от этого предупреждения, вы можете инициализировать fourier следующим образом:

fourier = np.zeros(k, dtype=complex)

Также формула для Дискретного преобразования Фурье включает суммирование по частотам, охватывающим полный [0,1) диапазон. Чтобы получить ДПФ с 1000 точками (как в вашем коде), вам нужно будет использовать

freq = np.arange(0,1,.001)

Это приведет к спектру, который включает 2 пика: один на ожидаемой частоте и другой симметричный выше частоты Найквиста. Обычно при построении спектра сигналов с действительным знаком отбрасывают результаты выше частоты Найквиста (но используйте полный спектр в своей функции IFT).

Наконец, как отметил GrimTrigger: