【问题标题】:Amplitude and phase spectrum. Shifting the phase leaving amplitude untouched幅度和相位谱。移动相位,保持幅度不变
【发布时间】:2018-09-05 07:53:08
【问题描述】:

我有在相关点具有等效间隔和相应测量值的数据。例如,这里是我拥有的数据的摘录:

y =[2.118, 2.1289, 2.1374, 2.1458, 2.1542, 2.1615, 2.1627, 2.165 2.1687...]

点之间的间隔是0.1

所以,我需要从数据中得到幅度谱(幅度与频率)以及相位谱(相位角与频率)。 此外,我应该将数据的相位移动负 90 度(-pi/2)。

在移相并保持幅度不变时,我需要进行反 fft 并获得新信号。我想在 Python 中做到这一点。

您能给我一个执行此操作的示例吗?

我使用的代码取自另一个 SO question,但我做了一些修改

## Perform FFT WITH SCIPY
signalFFT = np.fft.fft(y)

## Get Power Spectral Density
signalPSD = np.abs(signalFFT) ** 2
signalPhase = np.angle(signalFFT)

## Shift the phase py +90 degrees
new_signalPhase =(180/np.pi)*np.angle(signalFFT)+90 


## Get frequencies corresponding to signal 
fftFreq = np.fft.fftfreq(len(signalPSD), 0.1)

## Get positive half of frequencies
i = fftFreq>0

##
plt.figurefigsize=(8,4)
#plt.plot(fftFreq[i], 10*np.log10(signalPSD[i]));

plt.plot(fftFreq[i], new_signalPhase[i]);
plt.ylim(-200, 200);
plt.xlabel('Frequency Hz');
plt.ylabel('Phase Angle')
plt.grid()
plt.show()

问题是我想重新生成具有相同幅度但相移的信号。我知道答案与 ifft 有关,但我应该如何为它准备数据?您能否建议我进一步的步骤。

output

【问题讨论】:

  • 您好,Kamran,欢迎来到 SO!在这里提出问题时,重要的是向我们展示您已经尝试解决的问题,并非常具体地说明您遇到的问题。阅读how do I ask a good question?。这将帮助我们帮助您。目前,我们很难判断您是不了解 FFT 的工作原理,还是根本不知道如何编写代码。
  • 编写一些东西...一次只问一个问题
  • Kamran,请在下面查看我的答案和代码。相移信号由三行代码产生,傅里叶变换、乘以 exp(i 相位)和逆变换。
  • @DrM 非常感谢。如果我错了,请纠正我,乘以 exp(i phase) 只是简单地添加 exp(i w t) 的幂,我的意思是在 DFT 表达式中。我说对了吗?
  • exp( i w t ) x exp( i phase ) 等于 exp[ i (w t + phi ) ],即它改变了相位。换句话说,在将傅里叶变换信号乘以 exp( i 相位 ) 之后,您所拥有的是移相相同相位的信号的 FT。这适用于任何输入信号。

标签: python signal-processing fft ifft phase


【解决方案1】:

这是您的代码,经过一些修改。我们应用傅里叶变换,对变换后的信号进行相移,然后进行傅里叶逆变换以产生相移的时域信号。

请注意,变换是使用 rfft()irfft() 完成的,相移是通过简单地将变换后的数据乘以 cmath 来完成的.rect(1.,phase)。相移相当于将复数变换信号乘以exp(i * phase)。

在图形的左侧面板中,我们显示了原始信号和偏移信号。新信号提前 90 度。在右侧面板中,我们在左轴上显示功率谱。在这个例子中,我们有一个单一频率的功率。相位绘制在右轴上。但同样,由于我们只有一个频率的功率,相位谱在其他任何地方都会显示噪声。

#!/usr/bin/python

import matplotlib.pyplot as plt
import numpy as np
import cmath

# Generate a model signal
t0 = 1250.0
dt = 0.152
freq = (1./dt)/128

t = np.linspace( t0, t0+1024*dt, 1024, endpoint=False )
signal = np.sin( t*(2*np.pi)*freq )

## Fourier transform of real valued signal
signalFFT = np.fft.rfft(signal)

## Get Power Spectral Density
signalPSD = np.abs(signalFFT) ** 2
signalPSD /= len(signalFFT)**2

## Get Phase
signalPhase = np.angle(signalFFT)

## Phase Shift the signal +90 degrees
newSignalFFT = signalFFT * cmath.rect( 1., np.pi/2 )

## Reverse Fourier transform
newSignal = np.fft.irfft(newSignalFFT)

## Uncomment this line to restore the original baseline
# newSignal += signalFFT[0].real/len(signal)


# And now, the graphics -------------------

## Get frequencies corresponding to signal 
fftFreq = np.fft.rfftfreq(len(signal), dt)

plt.figure( figsize=(10, 4) )

ax1 = plt.subplot( 1, 2, 1 )
ax1.plot( t, signal, label='signal')
ax1.plot( t, newSignal, label='new signal')
ax1.set_ylabel( 'Signal' )
ax1.set_xlabel( 'time' )
ax1.legend()

ax2 = plt.subplot( 1, 2, 2 )
ax2.plot( fftFreq, signalPSD )
ax2.set_ylabel( 'Power' )
ax2.set_xlabel( 'frequency' )

ax2b = ax2.twinx()
ax2b.plot( fftFreq, signalPhase, alpha=0.25, color='r' )
ax2b.set_ylabel( 'Phase', color='r' )


plt.tight_layout()

plt.show()

这是图形输出。

【讨论】:

    猜你喜欢
    • 2017-11-13
    • 1970-01-01
    • 1970-01-01
    • 2021-05-25
    • 1970-01-01
    • 2011-09-17
    • 2011-12-07
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多