Амплитудный и фазовый спектр. Сдвиг фазы, оставляя амплитуду нетронутой

Question

У меня есть данные, которые имеют эквивалентные интервалы и соответствующие измерения в соответствующих точках. Например, вот выдержка из данных, которые у меня есть:

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, но я сделал некоторые изменения

## 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, но как мне подготовить данные для него? Не могли бы вы посоветовать мне дальнейшие шаги?

выход

Answer 1

Вот ваш код с некоторыми изменениями. Мы применяем преобразование Фурье, сдвигаем фазу преобразованного сигнала, а затем выполняем обратное преобразование Фурье для получения сдвинутого по фазе сигнала во временной области.

Обратите внимание, что преобразования выполняются с помощью rfft () и irfft () , и что фазовый сдвиг выполняется простым умножением преобразованных данных на cmath.rect (1., фаза) . Сдвиг фазы эквивалентен умножению комплексного преобразованного сигнала на 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()

А вот графический вывод.