How can I transform data in Time-Domain from excel to Frequency-Domain by using FFT in order to get PSD?

Hi everyone ,

I have data from exel in Time-Domain ( Time (ns) and real values) and i want to transform it to Frequency-domain by using FFT in order to get PSD . I have written below script and it doesn't work .

% data i got them Excel file , frist column is time(ns) and 2nd column real values. 
% length of data is 4760 .
time = data(:,1);       % sampling time 
signal = data(:,2);     % signal data in Time-Domain 
L=length(signal);       % Length of signal 
Ts=time, 
Fs=1/Ts ,               % sampling frequency 

NFFT = 2^nextpow2(L);   % Next power of 2 from length of y
Y = fft(signal,NFFT)/L;

f = Fs/2*linspace(0,1,NFFT/2+1)

figure(1),
plot(f,2*abs(Y(1:NFFT/2+1)));  % Plot single-sided amplitude spectrum.
grid on
title('Single-Sided Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
Y1 = fft(signal,NFFT)*Ts;
Y1 = fftshift(Y);
figure(2);
plot(f,abs(Y1));
grid on
title('Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y1(f)|');

figure(3);
plot(f,angle(Y));
grid on
title('Phase of y(t)');
xlabel('Frequency (Hz)');
ylabel('Phase of Y(f)');

Many thanks in advance for your time.

Answer 1

[Wayne King]

One mistake you're making is:

time = data(:,1);       % sampling time 
Ts = time;
Fs = 1/Ts;

The above does not give you the sampling frequency.

This will give you the proper sampling frequency

time = data(:,1);
Ts = mean(diff(time));
Fs = 1/Ts;

Then another mistake you make is when you go back and try to plot a two-sided Fourier transform against a frequency vector you constructed for just 1/2 the frequency interval. You can fix that with:

freq = -Fs/2+Fs/NFFT:Fs/NFFT:Fs/2;
 Ycent = fftshift(Y);
 plot(freq,abs(Ycent));

You've already obtained the Fourier transform once in your code, you don't need to do it again.

Then

   plot(freq,angle(Ycent))

I'm not going to quibble about your scaling for your PSD estimate, which does not technically yield a PSD estimate, but it's most likely fine for what you're doing.

If you have the Signal Processing Toolbox, why not use periodogram()? That will handle the proper scaling for you.

Answer 2

[gjhygf]

Fs = 1000;

t = 0:1/Fs:1-1/Fs;

signal = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

N = length(signal);

signal_fft = fft(signal);

P2 = abs(signal_fft/N);

P1 = P2(1:N/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(N/2))/N;

figure;

plot(f,P1);

title('Single-Sided Amplitude Spectrum of signal(t)');

xlabel('f (Hz)');

ylabel('|P1(f)|');

reconstructed_signal = ifft(signal_fft);

figure;

plot(t, signal);

hold on;

plot(t, real(reconstructed_signal), 'r');

legend('Original Signal', 'Reconstructed Signal');

xlabel('Time (s)');

ylabel('Amplitude');

title('Signal and its Reconstruction');

Answer 3

[gjhygf]

Fs = 1000;

t = 0:1/Fs:1-1/Fs;

signal = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

N = length(signal);

signal_fft = fft(signal);

P2 = abs(signal_fft/N);

P1 = P2(1:N/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(N/2))/N;

figure;

plot(f,P1);

title('Single-Sided Amplitude Spectrum of signal(t)');

xlabel('f (Hz)');

ylabel('|P1(f)|');

reconstructed_signal = ifft(signal_fft);

figure;

plot(t, signal);

hold on;

plot(t, real(reconstructed_signal), 'r');

legend('Original Signal', 'Reconstructed Signal');

xlabel('Time (s)');

ylabel('Amplitude');

title('Signal and its Reconstruction');