FFT of a signal based on simulation data

Hi all,

I'm a matlab newbie, and have been searching for an answer to this for a few days with no success. I've got some data from another program (CST microwave) and I need to FFT it. Here's the portion of my code that should be giving me a result:

fs = length(tr)/tr(end);
N = 1024;
Y = fft(xr,N);
k = [-N/2:N/2-1];
plot(k*fs/N,abs(fftshift(Y)))

However, the result I get seems wrong. Are there any obvious mistakes in this code? I'd be grateful for any insight anyone might have. Fs is the sampling frequency, and I'm trying to get it by using the simulation data's time vector (i.e: the signal I need to FFT lasts 40ns, and there are 4001 points).

Thanks for your time.

Answer 1

[Dr. Seis]

I usually do:

df = fs / N;
Nyq = fs / 2;
f = -Nyq : df : Nyq-df;
plot(f , abs(fftshift(Y)));

But that looks like the same thing you are doing. If you apply the same thing to an odd number of traces (i.e., N = 4001), then that flow will have to change a bit.

f = -Nyq : df : Nyq; % for odd N

Also, if tr(1) = 0, then:

fs = (length(tr)-1)/tr(end);

And if tr(1) ~= 0, then:

fs = (length(tr)-1)/(tr(end)-tr(1));

Answer 2

[David]

Thank you! I hadn't thought about the tr(1) = 0 bit.

The rest seems to lead to the same result, though. It could be that the answer I'm getting is correct, but it seems strange to me. If it's any help, here is my full code (amended with your suggestions):

X0 = dlmread('X0_test.txt');
t0 = X0(:,1);
x0 = X0(:,2); 

Tx = dlmread('Tx_test.txt');
t1 = Tx(:,1);
x1 = Tx(:,2);  

%sampling the probe signal
ti = [];
j=1;
for i=1:0.01:40;
ti(j) = i-1;
j=j+1;
end
xi = interp1(t0,x0,ti);  

%sampling the transmitted pulse
tm = [];
n=1;
for m=1:0.01:40;
tm(n) = m-1;
n=n+1;
end
xm = interp1(t1,x1,tm);  

%obtaining the reflected signal by subtracting the transmitted pulse from the probe signal
tr = ti;
xr = xi - xm;
figure;
plot(tr,xr);
grid;
xlabel('Time (ns)');
ylabel('Intensity (V/m)');
title('Reflected Pulse (from Re-sampled data)');  

%FFT of xr
figure;
grid;
fs = (length(tr)-1)/tr(end);
N = 1024;
Y = fft(xr,N);
k = [-N/2:N/2-1];
plot(k*fs/N,abs(fftshift(Y)))

The FFT result looks like this: http://s7.postimage.org/khqig5wmj/FFT.jpg

I suppose it could be that the FFT is being applied correctly, and the axes are not formatted properly. Anyway, thank you for your answer, and if you have any more insight it would be very helpful.

[David]

I forgot to mention: the signal originally contains frequencies from 0 to 4GHz.

Answer 3

[Dr. Seis]

A couple of helpful tips:

ti = 0 : 0.01 : 39;

Will do the same thing as:

ti = [];
j=1;
for i=1:0.01:40;
ti(j) = i-1;
j=j+1;
end

In order to get a closer look at the frequencies that contribute the most to your timeseries, as well as taking into consideration that the magnitude of your amplitude spectrum is symmetric, you could follow your "plot" with:

xlim([0 10]);

To only look at that portion.

Also, I should mention that your amplitudes are off by a factor of 1/fs. Your "plot" function should be:

plot(k*fs/N,abs(fftshift(Y))/fs);

This is because the FFT assumes your sample rate is unity (1 sample per second) when it is calculating the discrete integral. Let me know if you want more information about why this is. It is necessary if you plan to use the actual value of the amplitude for further calculations, but if you plan on only looking at the shape of your frequency spectrum it is not necessary.

[David]

Thank you for the tips, I'll add the factor of 1/fs, as I do need the amplitudes.