Реализация numpy .fft.fft () в Python - PullRequest
Новая
       31

Реализация numpy .fft.fft () в Python

6 голосов
/ 19 апреля 2020

Я пытаюсь получить спектр периодических сигналов c, используя функцию fft. Затем нанесите на график величину и фазу преобразования. Графики величины в порядке, но фазовые графики совершенно неожиданны. Например, я использовал функции Sin³ (t) и Cos³ (t). Код, который я использовал: 

import matplotlib.pyplot as plt
import numpy.fft as nf
import math
import numpy as np
pi = math.pi

N=512

# Sin³(x)
t=np.linspace(-4*pi,4*pi,N+1);t=t[:-1]
y=(np.sin(t))**3
Y=nf.fftshift(nf.fft(y))/N
w=np.linspace(-64,64,N+1);w=w[:-1]

plt.figure("0")
plt.subplot(2,1,1)
plt.plot(w,abs(Y),'ro',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"$|Y|$",size=16)
plt.title("Spectrum of sin\u00B3(t)")
plt.grid(True)
plt.subplot(2,1,2)
ii=np.where(abs(Y)>1e-3)
plt.plot(w[ii],np.angle(Y[ii]),'go',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"Phase of $Y$",size=16)
plt.xlabel(r"$\omega$",size=16)
plt.grid(True)

# Cos³(x)

t=np.linspace(-4*pi,4*pi,N+1);t=t[:-1]
y=(np.cos(t))**3
Y=nf.fftshift(nf.fft(y))/N
w=np.linspace(-64,64,N+1);w=w[:-1]

plt.figure("1")
plt.subplot(2,1,1)
plt.plot(w,abs(Y),'ro',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"$|Y|$",size=16)
plt.title("Spectrum of cos\u00B3(t)")
plt.grid(True)
plt.subplot(2,1,2)
ii=np.where(abs(Y)>1e-3)
plt.plot(w[ii],(np.angle(Y[ii])),'go',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"Phase of $Y$",size=16)
plt.xlabel(r"$\omega$",size=16)
plt.grid(True)

Полученные графики:

i) Для Sin³ (t) - Величина и фаза спектра Sin³ (t)

ii) Для Cos³ (t) - Величина и фаза спектра Cos³ (t)

Как видно из приведенных выше ссылок, величина обоих функции в порядке. Фаза Спектра Sin³ (t) является правильной, как и ожидалось.

Поскольку Cos³ (t) является действительным и четным, а расширение в терминах сложных экспонент показывает, что Фаза Коэффициентов равна 0. Но график График показывает совершенно другой ответ (см. 2 ): при w = -3 фаза составляет около 5 радиан. Какую ошибку я совершил и как правильно реализовать БПФ.

1 Ответ

3 голосов
/ 19 апреля 2020

np.fft ведет себя как ожидалось; это ваш заговор, который вызывает путаницу. Звонок на plt.tight_layout() должен помочь вам разобраться:

Spectrum of cos t ^ 3

Если вы внимательно посмотрите на ось y фазы для cos³ (t) вы увидите, что все значения имеют коэффициент 1e-17. Таким образом, при w = -3 фаза равна , а не около 5 радиан, фактически она равна 5e-17 радианам (что для всех намерений и целей равно нулю).

Я убрал ваш код, пока Я смотрел в это: 

import matplotlib.pyplot as plt
import numpy.fft as nf
import numpy as np

plt.rcParams['axes.grid'] = True

pi = np.pi
N = 512

t = np.linspace(-4*pi, 4*pi, N+1)[:-1]


def fft_func(x, func, N):
    y = func(x) ** 3
    Y = nf.fftshift(nf.fft(y)) / N
    w = np.linspace(-64, 64, N+1)[:-1]
    return y, Y, w


def plot_mag_phase(w, Y, func_name):
    fig, ax = plt.subplots(2, 1)

    ii = np.where(abs(Y) > 1e-3)

    ax[0].plot(w, abs(Y), 'ro', lw=2)
    ax[1].plot(w[ii], np.angle(Y[ii]), 'go', lw=2)

    ax[0].set_xlim(-4, 4)
    ax[1].set_xlim(-4, 4)

    ax[0].set_ylabel("$|Y|$", size=16)
    ax[1].set_ylabel("Phase of $Y$", size=16)

    ax[0].set_title("Spectrum of {}\u00B3(t)".format(func_name))
    ax[1].set_xlabel(r'$\omega$', size=16)

    fig.tight_layout()

# sin x ^ 3
y, Y, w = fft_func(t, np.sin, N)
plot_mag_phase(w, Y, 'sin')

# cos x ^ 3
y, Y, w = fft_func(t, np.cos, N)
plot_mag_phase(w, Y, 'cos')
...