After the FFT(Frequency-Amplitude) of the seismic wave(Time-Acceleration), I want to reconstruct the seismic wave(Time-Acceleration) by IFFT again.
5 次查看(过去 30 天)
显示 更早的评论
First, FFT was performed using time and acceleration data.
And reconstruct the seismic wave(Time-Acceleration) by IFFT again.
%Fast Fourier Spectrum
Time=[0, 0.01, 0.02, 0.03, ..., 100 ] % Time (s)
Acceleration=[0.0001, 0.0005, 0.0006, ..., 0.6842 ] % Acceleration (m/s^2)
figure;
plot(Time,Acceleration);
grid;
xlabel('Time (sec)');
ylabel('Acceleration(m/s^2)');
title('Input Ground Motions Time History');
n = length(Acceleration);
Ts = Time(2)-Time(1); % Sample Time
Fs = 1/Ts; % Sampling Frequency
NFFT = 2^nextpow2(n); % Next power of 2 from length of data
Y = fft(Acceleration,NFFT)/n;
f = Fs/2*linspace(0,1,NFFT/2+1);
Iv = 1:length(f); % Index Vector
%Filtering and Combine
combined=Y(Iv); % Fitering and Combined - Multiply the Fourier amplitude by a certain expression.
% For example Y(Iv).*exp(0.011.*f),
% But Now, We do not multiply the expression because we want to output the input value as it is.
%IFFT
COM=ifft(combined,NFFT)/n;
figure;
plot(COMreal);
grid;
xlabel('Time (sec)');
ylabel('Acceleration(m/s^2)');
title('IFFT Ground Motions Time History');
There is a problem here.
After performing the FFT, even though no formula was multiplied, the input time history and the IFFT time history are different when the IFFT is performed.
It can be seen that the peak value is different(input=4.04, output=2.02*10^-8) and the x-axis (Time) value is longer.
2 个评论
Bjorn Gustavsson
2020-5-7
When you plot your final panel you do
plot(COMreal);
Assuming that COMreal is the real part of COM you don't plot it as a function of time only against the array index. It seems to line up OK to me, try:
plot(Time,real(COM));
回答(2 个)
Mehmed Saad
2020-5-7
编辑:Mehmed Saad
2020-5-7
COM=ifft(Y,NFFT)*n;% divide in fft, multiply in ifft
figure;
COMreal = COM(1:length(Time)); % pick the data equal to actual data size
plot(Time,COMreal)% plot against time if you just write plot(COMreal)
% it will be plotted against samples
grid;
xlabel('Time (sec)');
ylabel('Acceleration(m/s^2)');
title('IFFT Ground Motions Time History');
0 个评论
David Goodmanson
2020-5-7
Hello haeyeon,
A couple of things here. First of all, you are using nextpow2 to go from what looks like 10001 points to 16384. That is going to cause some issues with the time array in the ifft calculation and with the amplitude. You can adjust for that, but it is much easier to just not bother with nextpow2 and go with the original number of points. The fft will be plenty fast anyway. So suppose that n = 10001 or whatever it actually is.
To go to the frequency domain, you are dividing the fft by n. I'm not sure if that is the standard way of doing things in seismic calculatons, but I assume here that it is. That means that when you go back to the time domain, you have to multiply the ifft by n in order to make up for dividing by n. You are dividing the ifft by n, which means it's off by a factor of n^2 which is 10^8. Bringing the second plot up by that factor brings the agreement to a factor of 2.
For the time axis in the second plot you are using not the time, but just the index of the array, 1 to 16384. One advantage of not messing things up with nextpow2 is that for the second plot you can just use the original time array that you started with since the number of points has not changed.
That leaves the factor of 2 difference in amplitude, which is most likely due to truncating the frequency array to NFFT/2+1 points and then padding it out again with iffft(...,NFFT).
If you want to get back what you started with, the easiest way is
n = length(Acceleration);
Ts = Time(2)-Time(1); % Sample Time
Fs = 1/Ts; % Sampling Frequency
Y = fft(Acceleration)/n;
%IFFT
COM=ifft(Y)*n;
figure;
plot(t,real(COM)); % use real to eliminate small nuisance imaginary part
grid;
xlabel('Time (sec)');
ylabel('Acceleration(m/s^2)');
title('IFFT Ground Motions Time History');
this keeps all of Y so the ifft is going to work.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Earthquake Engineering 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!