Don´t know how to perform IFFT considering my frequency.
4 views (last 30 days)
Show older comments
Andriy Voshchenko
on 19 Sep 2022
Commented: Andriy Voshchenko
on 21 Sep 2022
Hi! I have a frequency domain data which I want to transform into time domain. I have the following file where the first column is the frequency in range [0 Nyquist frequency], second column is real component of the response and the 3rd colums is the imaginary part.
I know I have to use "ifft" function but I just don´t understand how to define it properly. I know the result should be time response doman with maximum value around 600-800 but when I launch the script I obtain peak values of 30. As far as I see there should be some way to input my frequency as the parameter. The help section includes this:
X = ifft(Y,n) returns the n-point inverse Fourier transform of Y by padding Y with trailing zeros to length n.
But I don´t understand what is the correlation between parameter "n" and my frequency, how do I know how many points to solve?
Part of my code as well as the attached data, keep in mind I have one sided spectra:
Mreaction=readmatrix("Mreaction_damped0.xlsx");
freq=Mreaction(:,1); % frequency of my response spectrum
hMoment=Mreaction(:,2)+i*Mreaction(:,3); % obtain complex number for the response moment
myTdata=ifft(hMoment,"symmetric"); % I´ve also tried with 1000 (and some others) points but the results just don´t
% make any sense to me
plot(myTdata) % Plot without time vector since I struggle to define it
Thanks in advance!
0 Comments
Accepted Answer
Matt J
on 19 Sep 2022
Edited: Matt J
on 19 Sep 2022
But I don´t understand what is the correlation between parameter "n" and my frequency, how do I know how many points to solve?
The key relationship involving n is,
n=1/(dF*dT)
where dT and dF are the time and frequency sampling intervals . The choice of n depends on what dF you have (in your case dF=0.03) and what dT you want. Also, n must be an integer.
When using fft(x,n) you are telling the code to append zeros to x making it length n. However, in applications to time-frequency analysis (like yours) I find it is often better just to manually pad x to the length that you want, because the padding has to be done in a very specific way in conjunction with fftshift.
I illustrate all this below:
Mreaction=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1129380/Mreaction_damped0.xlsx');
dT0=0.01; %maximum desired time sampling interval
dF=Mreaction(2)-Mreaction(1); %given frequency sampling interval
hMoment=complex( Mreaction(:,2) , Mreaction(:,3) );
n0=numel(hMoment);
n=max( ceil(1/dT0/dF) , n0 ); %ensures dT<=dT0=0.01
dT=1/dF/n; %actual time sampling interval
nAxis=( (0:n-1)-ceil((n-1)/2) ); %normalized axis
freqAxis=dF*nAxis;%frequency axis
timeAxis=dT*nAxis;%time axis
leftside=flip(conj(hMoment));
hMoment=[leftside(1:end-1); hMoment]; %two sided spectrum
hMoment(end+1:n)=0;
hMoment=fftshift(circshift(hMoment,1-n0)); %zero padding circulantly is complicated
myTdata=fftshift(real( ifft( ifftshift(hMoment)) ))/dT;
figure
plot(freqAxis,abs(hMoment));
xlabel Frequency
figure
plot(timeAxis, myTdata)
xlabel Time
5 Comments
Matt J
on 20 Sep 2022
Edited: Matt J
on 20 Sep 2022
There's no way to account for the scaling fully without seeing the code you used to to obtain the spectrum.
My code assumes that the given spectrum was correctly scaled to approximate the continuous Fourier transform. You would have needed to multiply the FFT by dT in order to ensure this.
Assuming you did so, then dividing the IFFT by dT is the appropriate way to invert. You can see in my code that I have done this, whereas Paul has not.
If you didn't scale the forward FFT, then Paul's scaling is probably correct and your expectations are probably wrong.
More Answers (1)
Paul
on 20 Sep 2022
Hi Andriy
I would attack the problem this way, assuming throughout the original time-domain signal is real and the freq vector in the file is in Hz.
Read in the data
%Mreaction=readmatrix("Mreaction_damped0.xlsx");
Mreaction = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1129380/Mreaction_damped0.xlsx');
freq = Mreaction(:,1); % frequency of my response spectrum
hMoment = Mreaction(:,2) + 1i*Mreaction(:,3); % obtain complex number for the response moment
Check the first point in the freq vector
freq(1)
It's not zero like it should be, so I'll synthesize a point at f = 0.
freq = [0;freq];
hMoment = [real(hMoment(1));hMoment];
The last point in the one-sided spectrum has a signficanct imaginary part
hMoment(end)
If this point truly corresponded to the Nyquist frequency, then it would be real and N, the length of the time-domain signal, would be even. But it's not, so it looks like N must be odd.
N = 2*numel(freq) + 1;
The indices for the signal and the difference between successive frequency points are
n = 0:(N-1);
dF = 0.03; % assuming this is Hz?
Therefore, the sampling frequency and sampling period are
Fs = N*dF;
Ts = 1/Fs;
The time domain signal is then
myTdata = ifft(hMoment,N,"symmetric");
plot(n*Ts,myTdata)
xlabel('Time (sec)')
I know that you were expecting the signal to have a maximum value in the 600-800 range. The appropriate scaling, if any, would depend on how the original frequency domain data was collected and scaled in the first place.
3 Comments
Paul
on 21 Sep 2022
As already stated, the appropriate scaling of the output of ifft() depends very much (entirely?) on how the frequency domain data was developed in the first place. Was it based on an FFT algorithm? If so, what scaling was assumed on that FFT? 1? 1/N? As @Matt J hypothesized, maybe the data was scaled to be consistent with the Continuous Time Fourier Transform, in which as a divide by Ts is in order as he showed, and seems to yield a result consistent with your expectations.
I agree with @Matt J that there really is not any way to know for sure without knowing exactly the process used to generate the frequency domaing data in the .xlsx file.
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!