11746

Scipy FFT - how to get phase angle

<h3>Question</h3>

I'm having trouble getting the phase of a simple sine curve using the scipy fft module in python. I followed this tutorial closely and converted the matlab code to python. However, no matter what phase I use for the input, the graph always shows 3. What am I missing?

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()

EDIT: Here is some simpler code. Still can't seem to get the phase.

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)]

The frequencies I get are clustered around 2*np.pi/7. But the phases I get are:

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

While there should be no phase at all.


<h3>Answer1:</h3>

This is the simplest code to show how to get the angles.

Note that I've created the signal y such that there is an integer number of periods in it (as I suggested in a comment, and @hotpaw2 suggested in their answer).

This is the problem with OP's code.

I used linspace to create the time axis, using the endpoint=False option. This is important, if t were to contain the number 10, then we no longer have an exact integer number of periods. When working with the discrete Fourier transform, it helps to think of what happens when you repeat the signal. Simply take y, and concatenate a copy of itself: np.hstack((y,y)). Is this new signal still a sampling of a single sine wave, or did we create something more complex? What happens at that point where the two copies meet?

<pre class="lang-py prettyprint-override">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()


<h3>Answer2:</h3>

An FFT measures circular phase, referenced to both the very beginning and very end of the input data window. If your input sine wave isn't exactly integer periodic in the FFT aperture, then there will be a discontinuity between the phase at the beginning and end of the window, thus the FFT phase measurement won't be what you might expect.

Instead, I recommend one begin (or reference) the input sinewave at the desired phase in the middle of your data window (at N/2), not the beginning. Then do an circular fftshift before the FFT. The resulting FFT phase measurement will then represent phase with respect to the middle of your original data window, where there is no discontinuity; and atan2() will represent the ratio between the even and odd components of your continuous input waveform as is more typically expected.

来源:https://stackoverflow.com/questions/54454723/scipy-fft-how-to-get-phase-angle

Recommend

  • VS Code extension: Hide commands from command palette
  • Can async-await be available in other .NET languages besides C#?
  • Get VS2008 to “tree-indent” partial classes (like code-behind files)
  • Unable to count the number of matches in Vim
  • FB JS-SDK cannot detect user logging out of FB directly
  • Is there a Python library to list primes?
  • Sum the content of 2 list in Groovy
  • Alternate locations for javax.comm.properties in Windows
  • What do folks use app/services/ in rails applications
  • having a condition to fragment adding
  • Dart HttpClient with Basic Authentication issues after flutter upgrade
  • Python multiprocessing: can I reuse processes (already parallelized functions) with updated global v
  • org.springframework.core.convert.ConverterNotFoundException: No converter found capable of convertin
  • How to add a progress ring to the splash screen in Windows 8?
  • How to set title name of the pdf. While viewing the Document(New Tab)
  • React / Webpack - “Module parse failed: Unexpected token - You may need an appropriate loader to han
  • web2py: How to execute instructions before delete using SQLFORM.smartgrid
  • get all files in git diff in intellij
  • How to use CoreFoundation in QuickTime SDK for Windows?
  • Vue.js 2: Vue cannot find files from /assets folder (v-for)
  • date changes on export kendoGrid
  • Refresh JSF component after custom javascript Ajax call
  • Complex multiple if statements
  • Slick: How can I combine a SQL LIKE statement with a SQL IN statement
  • Issues with converting data into a matrix after running lapply()
  • In metro, get all inherited classes of an (abstract) class?
  • About global variables in Node.js
  • 'float' object cannot be interpreted as an integer
  • How to change the host IP sent in emails to new GitLab users to a publicly visible IP, not the local
  • How to use Typescript with libraries like Ampersand.js that parse configs to build prototypes
  • How to use array in autohotkey?
  • Wireshark Display Filter for Unique Source/Destination IP and Protocol
  • Ember.js + JQuery-UI Tooltip - Tooltip does not reflect the model / controller changes
  • Year over Year Stats from a Crossfilter Dataset
  • Can I read another applications memory?
  • How can I ssh into a server that requires 2 password authentication using python's paramiko mod
  • Cloud Code: Creating a Parse.File from URL
  • How to integrate angular2-material (alpha 8.2) with angular2-Quickstart app
  • How to get rgb from transparent pixel in js