У меня есть данные с эквивалентными интервалами и соответствующими измерениями в соответствующих точках. В качестве примера приведу отрывок из имеющихся у меня данных:
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()
Проблема в том, что я хочу восстановить сигнал с теми же амплитудами, но со сдвигом фазы. Я знаю, что ответ связан с чем-то, связанным с ошибкой, но как мне подготовить для этого данные? Не могли бы вы посоветовать мне дальнейшие шаги.
Придумайте что-нибудь ... Задавайте только один вопрос за раз
Кямран, пожалуйста, посмотрите мой ответ и код ниже. Сигнал с фазовым сдвигом создается тремя строками кода, преобразованием Фурье, умножением на exp (фаза i) и обратным преобразованием.
@DrM большое спасибо большое. Поправьте меня, если я ошибаюсь, умножение на exp (i-фаза) - это просто сложение в степени exp (i w t), я имею в виду в выражении DFT. Я правильно понимаю?
exp (i w t) x exp (i phase) равно exp [i (w t + phi)], то есть сдвигает фазу. Другими словами, после умножения сигнала, преобразованного Фурье, на exp (i-фаза), у вас есть FT сигнала, сдвинутого на ту же самую фазу. Это работает для любого входного сигнала.
Попробуйте, и вы убедитесь, что это работает. И, пожалуйста, не забудьте отметить это правильно и / или проголосовать за него.
@DrM есть ли способ не нормализовать newSignal? Кажется, что новостной сигнал помещается около 0 и не соответствует исходным значениям. Я приложил вывод в описании проблемы выше.
Если входной сигнал имеет смещение по постоянному току, оно должно увеличиваться в элементе нулевой частоты. Я отредактировал код в своем ответе, чтобы показать, как добавить это обратно, чтобы восстановить исходное смещение. Раскомментируйте новую строку, чтобы увидеть, как она работает.
Вот ваш код с некоторыми изменениями. Мы применяем преобразование Фурье, сдвигаем по фазе преобразованный сигнал, а затем выполняем обратное преобразование Фурье, чтобы получить сдвинутый по фазе сигнал во временной области.
Обратите внимание, что преобразования выполняются с помощью 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()
А вот и графический вывод.
Привет, Кямран, и добро пожаловать в ТАК! Задавая здесь вопрос, важно показать нам, что вы уже пытались решить для своих проблем, и очень конкретно указать, с чем вы застряли. Прочтите как мне задать хороший вопрос?. Это поможет нам помочь вам. Сейчас нам трудно сказать, не понимаете ли вы, как работает БПФ, или просто не знаете, как его кодировать.