Amplitude and phase spectrum. Shifting the phase leaving amplitude untouched

Kamran Abbasov picture Kamran Abbasov · Sep 5, 2018 · Viewed 7.3k times · Source

I have the data that has equivalent intervals and corresponding measurements at relevant points. As an example, here is the excerpt of the data I have:

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

interval between the points is 0.1

So, what I need to get out of the data is the Amplitude spectrum (amplitude vs frequency) and also phase spectrum (phase angle vs frequency). In addition I should shift the phase of the data by negative 90 degrees (-pi/2).

Upon shifting the phase and leaving the amplitude untouched, I need to do the inverse fft and get the new signal. I want to do this in Python.

Could you please give me an example of performing this.

The code that I have used, was taken from another SO question, but I have done some modifications

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

The problem is that I want to regenerate the signal, with the same amplitudes but shifted phase. I know that the answer is smth related to ifft, but how should I prepare data for it? Could you please advice me the further steps.

output

Answer

DrM picture DrM · Sep 12, 2018

Here is your code with some modifications. We apply a Fourier transform, phase shift the transformed signal, and then perform the inverse Fourier transform to produce the phase shifted time domain signal.

Notice that the transforms are done with rfft() and irfft(), and that the phase shift is done by simply multiplying the transformed data by cmath.rect(1.,phase). The phase shift is equivalent to multiplying the complex transformed signal by exp( i * phase ).

In the graphics, in the left panel, we show the original and shifted signals. The new signal is advanced by 90 degrees. In the right panel we show the power spectrum on the left axis. In this example we have power at a single frequency. Phase is graphed on the right axis. But again, since we have power at only one frequency, the phase spectrum shows noise everywhere else.

#!/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()

And here is the graphical output.

Raw and shifted signals, power and phase spectrum