Scipy FFT - как получить фазовый угол - PullRequest
Новая
       28

Scipy FFT - как получить фазовый угол

0 голосов
/ 31 января 2019

У меня проблемы с получением фазы простой синусоиды с использованием модуля scipy fft в python.Я внимательно следил за этим учебником и преобразовал код matlab в python.Однако независимо от того, какую фазу я использую для ввода, график всегда показывает 3. Чего мне не хватает? 

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import cmath
A=10
fc = 10
phase=60
fs=32#Sampling frequency with oversampling factor of 32

t = np.arange(0,2,1/fs)

#Convert the phase shift to radians from degrees.
phi = phase*np.pi/180

x=A*np.cos(2*np.pi*fc*t+phi)

N=256
X = scipy.fftpack.fftshift(scipy.fftpack.fft(x,N))/N

df=fs/N #Frequency resolution.
sampleindex = np.arange(-N/2,N/2,1) #Ordered index for FFT plot.
f = sampleindex*df #x-axis index continued to ordered frequencies

raw_phases = np.angle(X)

X2=np.copy(X)#Store the FFT results in another array.
#Detect very small numbers and ignore them.
tau = max(abs(X))/10
X2[abs(X)<tau]=0

phase=[cmath.phase(i) for i in X2]
plt.plot(f,phase)
plt.show()

enter image description here

РЕДАКТИРОВАТЬ: Здеськакой-то более простой кодКажется, до сих пор не получается получить фазу. 

y = 28*np.sin(2*np.pi/7*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
phase = np.angle(yf)
yf = np.abs(yf[:N//2])
phase = phase[:N//2]
xf = xf[1:]
yf = yf[1:]
phase = phase[1:]
yf = yf-np.mean(yf)
#The frequencies seem to always be scaled by 0.1433, not sure why.
c = 2*np.pi/7/0.143301
freqs = xf[yf>np.std(yf)]*c
phases = phase[yf>np.std(yf)]

Частоты, которые я получаю, сгруппированы вокруг 2 * np.pi / 7.Но фазы, которые я получаю: 

array([-0.217795  , -0.22007488, -0.22226087,  2.91723935,  2.91524011,
    2.91333367])

Пока не должно быть никакой фазы.

Ответы [ 2 ]

0 голосов
/ 31 января 2019

Это самый простой код, показывающий, как получить углы.

Обратите внимание, что я создал сигнал y так, что в нем есть целое число периодов (как я предложил * 1004).* в комментарии , а @ hotpaw2 предложил в своем ответе ).

Это проблема с кодом OP.

Я использовал linspace для созданияось времени, используя опцию endpoint=False.Это важно, если t будет содержать число 10, то у нас больше не будет точного целого числа периодов.При работе с дискретным преобразованием Фурье помогает думать о том, что происходит, когда вы повторяете сигнал.Просто возьмите y и объедините его копию: np.hstack((y,y)).Этот новый сигнал все еще является выборкой синусоидальной волны, или мы создали что-то более сложное?Что происходит в тот момент, когда две копии встречаются? 

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

phase = np.pi / 4
t = np.linspace(0, 10, num=200, endpoint=False)
y = np.cos(2 * np.pi * t + phase)
Y = scipy.fftpack.fftshift(scipy.fftpack.fft(y))
f = scipy.fftpack.fftshift(scipy.fftpack.fftfreq(len(t)))

p = np.angle(Y)
p[np.abs(Y) < 1] = 0
plt.plot(f, p)
plt.show()

plot of the phase of the Frequency spectrum of a pure sine wave

0 голосов
/ 31 января 2019

БПФ измеряет круговую фазу, относящуюся как к самому началу, так и к самому концу окна входных данных.Если входная синусоида не является целочисленной периодической в ​​апертуре БПФ, то между фазой в начале и в конце окна будет разрыв, поэтому измерение фазы БПФ будет не таким, как вы могли ожидать.

Вместо этого я рекомендую начинать (или ссылаться) синусоиду ввода в нужной фазе в середине окна данных (в N / 2), а не в начале.Затем сделайте круговой сдвиг перед БПФ.Полученное в результате измерение фазы БПФ будет представлять фазу относительно середины исходного окна данных, где нет разрыва;и atan2 () будет представлять соотношение между четными и нечетными компонентами вашего непрерывного входного сигнала, как обычно ожидается.

...