applying ifft to a frequency domain data

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
ouess
User
Beiträge: 3
Registriert: Mittwoch 8. August 2018, 11:20

Hello,

I am trying to rebuild a time signal domain from frequency domain data.

For this I create a time domain data and then I apply the fft (this part is only to create a random signal and get frequency domain equivalent signal). Like this I have the input for my script which is supposed to rebuild a signal from frequency domain data.

The ifft is giving me an acceptable signal at the start of the time line but it starts to shift when the time is increasing. Could you please tell me what is my mistake ?

I am sorry I don't speak German.

Code: Alles auswählen

import numpy as np 
from scipy import pi 
import matplotlib.pyplot as plt
#%matplotlib inline
from scipy.fftpack import fft 

min=0
max=400

def calculateFFT (timeStep,micDataX,micDataY):
	
	n=micDataX.size
	
	FFT=np.fft.fft(micDataY)
	fft_amlitude=2*abs(FFT)/n
	fft_phase=np.angle(FFT)
	fft_freq= np.fft.fftfreq(n, d=timeStep) #not used created manually (7 lines) check pi_fFreqDomainCreateConstantBW it is kept here to compare sizes 
	upper_bound=int((n)/2)
	
	return fft_freq[1:upper_bound],fft_amlitude[1:upper_bound],fft_phase[1:upper_bound]

	
def calculateI_FFT (n,amplitude_spect,phase_spect):
	
	data=list()
	
	for mag,phase in zip(amplitude_spect,phase_spect):
		data.append((mag*n/2)*(np.cos(phase)+1j* np.sin(phase)))
	full_data=list(data)
	i_data=np.fft.irfft(data)
	return i_data

#sampling rate and time vector
start_time=0 #sec
end_time= 2 
sampling_rate=1000 #Hz
N=(end_time-start_time)*sampling_rate 

#Freq domain peaks
peak1_hz=60 # freq of peak 
peak1_mag= 25 

peak2_hz=270 # freq of peak 
peak2_mag= 2 

#Vibration data generation
time =np.linspace(start_time,end_time,N)
vib_data=peak1_mag*np.sin(2*pi*peak1_hz*time)+peak2_mag*np.sin(2*pi*peak2_hz*time)

#Data plotting 
plt.plot(time[min:max],vib_data[min:max])

# fft 
time_step=1/sampling_rate
fft_freq,fft_data,fft_phase=calculateFFT(time_step,time,vib_data)

#ifft 
i_data=calculateI_FFT(N,fft_data,fft_phase)
	
#plotting
plt.plot(time[min:max],i_data[min:max])
plt.xlabel("Time (s)")
plt.ylabel("Vibration (g)")
plt.title("Time domain")


plt.show()

Thank you
Wess
Benutzeravatar
ThomasL
User
Beiträge: 1366
Registriert: Montag 14. Mai 2018, 14:44
Wohnort: Kreis Unna NRW

Hi Wess

Code: Alles auswählen

#Vibration data generation
print(N) # N == 2000
time = np.linspace(start_time, end_time, N)
print(time.shape) # gives (2000,)
vib_data = peak1_mag * np.sin(2 * pi * peak1_hz * time) + peak2_mag * np.sin(2 * pi * peak2_hz * time)
print(vib_data.shape) # gives (2000,)

#Data plotting
plt.plot(time[minimum:maximum], vib_data[minimum:maximum])

# fft
time_step = 1 / sampling_rate
fft_freq, fft_data, fft_phase = calculateFFT(time_step, time, vib_data)

#ifft
i_data = calculateI_FFT(N, fft_data, fft_phase)
print(i_data.shape) # gives (1996,)
I would assume that i_data should also be of shape (2000,)?
Ich bin Pazifist und greife niemanden an, auch nicht mit Worten.
Für alle meine Code Beispiele gilt: "There is always a better way."
https://projecteuler.net/profile/Brotherluii.png
ouess
User
Beiträge: 3
Registriert: Mittwoch 8. August 2018, 11:20

Hi ThomasL,

Danke !
Yes you are right, there is something wrong with what I am doing.
For the moment I am almost sure it comes from the components A[0] and the negative frequency values (symetrical) but I can't find how to do it.
Benutzeravatar
ThomasL
User
Beiträge: 1366
Registriert: Montag 14. Mai 2018, 14:44
Wohnort: Kreis Unna NRW

Your function calculateFFT should return arrays starting at index 0 not 1

Code: Alles auswählen

def calculateFFT (timeStep, micDataX, micDataY):
    .........
    return fft_freq[0:upper_bound], fft_amlitude[0:upper_bound], fft_phase[0:upper_bound]
Are you possibly working with MATLAB also?
I know that in MATLAB arrays/matrices start at index 1
Ich bin Pazifist und greife niemanden an, auch nicht mit Worten.
Für alle meine Code Beispiele gilt: "There is always a better way."
https://projecteuler.net/profile/Brotherluii.png
ouess
User
Beiträge: 3
Registriert: Mittwoch 8. August 2018, 11:20

This is the Zero frequency term you won't find it in frequency dependent signal no ?? so I am skipping it for FFT but when you need to go back to time domain it looks like it is needed now I need to figure out if it is possible to calulate it or do without it.
Antworten