Trying to verify FFT code for research. Verification isn't working however.

6 views (last 30 days)
Hello!
I have been working on a Matlab project for some time, involving some seismic wave analysis using FFT in Matlab. The results aren't as we had hoped, so I have had to verify that the matlab FFT function does as expected of it. From what I can tell I use it the exact same way as the FFT Matlab guide (<http://www.mathworks.com/help/matlab/ref/fft.html)>, with its example being:
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
To test it, I am using equations 5,6,7 on this paper (<http://www.seg2.ethz.ch/dalguer/courses/earthquake/Beresnev2003BSSA.pdf>) comparing the FFT of the "triangle" function to the squared sinc function. Presumably, if it is done correctly, they should be equal if given the same parameters, but what I am finding is that they are not equal.
Attached below is the code to test these equations that I wrote up.
clear all
close all
U = 1; %Just a scalar.
dt = .01;
t0 = 3; %Because why not
t1 = 0:dt:t0/2; %For the first half.
t2 = t0/2+dt:dt:t0; %For the second half. Make sure to not have two points repeat.
%Just for curiosity, the other base equation. Does it do what it should?
firsth = 2*U*(t1/t0).^2;
sech = U*[1-2*(1-t2/t0).^2];
bothh = [firsth' ; sech']'; %put them togehter.
%plot(bothh) %Does it plot? Yes! Max of 1 too, when t=t0.
firsthalf = 4*(U/t0^2)*t1;
secondhalf = 4*(U/t0^2)*(t0-t2);
bothhalves = [firsthalf' ; secondhalf']';
%figure()
%plot(bothhalves) %Yep, definitely a nice triangle right there.
DATA2=detrend(bothhalves,'linear');
Fs=1/dt;
l=length(bothhalves);
L=l;
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(DATA2,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1); %Gives me the frequencies
y2 = abs(Y(1:NFFT/2+1)); %Gives me the absolute value of the magnitudes per frequency
omeega = 2*pi*f; %Just a definition.
checkit = U*[(sin(omeega*t0/4))./(omeega*t0/4) ].^2; %The numbers.
%The only reason the next part exists is because the first unit of the sine
%function gives me a NaN. So, when gettin the difference, it is also a NaN.
%The first point for y2 is basically 0, something*10^-17. So I got rid of
%the first point of both so that I could plot the difference between the
%two.
checkit(1) = [];
y2(1) = [];
f(1)=[];
dif = y2./checkit;
y2 = y2*dt;
figure()
loglog(f,checkit,'g') %sine function
hold on
loglog(f,y2,'r') %fft
figure
plot(dif)
Attached is the picture of the two, and clearly they are similar but there are great difference between the two. Early on the orders of magnitude are extremely far apart, which might be causing error in my actual project, as it is using the matlab FFT the same way.
Any insight as to why they are different? If not, any ideas of different tests for my code to insure accuracy? Maybe I am missing something pretty obvious. Thanks in advance!

Answers (1)

Honglei Chen
Honglei Chen on 11 Aug 2014
I couldn't find the pictures but my guess is very likely you are off by the magnitude of dt because you are comparing a discrete Fourier transform result with a continuous Fourier transform result.
  1 Comment
Nick
Nick on 14 Aug 2014
Sorry, first time posting. I believe I have it attached now to this response. I don't know if a result of just magnitude different would be a result of what you said.
And if it is, could you explain in a bit more detail as to why it is so? Thanks!

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!