FFT of a signal based on simulation data

7 次查看（过去 30 天）
David
David 2012-1-17
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.

请先登录，再回答此问题。

采纳的回答

Dr. Seis
Dr. Seis 2012-1-17
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));

更多回答（2 个）

David
David 2012-1-17
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.
  1 个评论
David
David 2012-1-17
I forgot to mention: the signal originally contains frequencies from 0 to 4GHz.

请先登录，再进行评论。

Dr. Seis
Dr. Seis 2012-1-17
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.
  1 个评论
David
David 2012-1-18
Thank you for the tips, I'll add the factor of 1/fs, as I do need the amplitudes.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Signal Operations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by