MATLAB: How to transform data in Time-Domain from excel to Frequency-Domain by using FFT in order to get PSD

fftmatlab coder

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.

Best Answer

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.

Related Solutions

MATLAB: I get the peaks at the wrong frequency and with the wrong amplitude.

hello Anna

look at my fft function at the end of my code

clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
data = importdata('SFC5S_nov25_ST_1_1.txt');
dt = 1e-6; % 1 micro seconds
signal = data(:,1);
samples = length(signal);
Fs = 1/dt; % sampling frequency (Hz)
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.95;
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed

if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end
figure(1),plot(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

MATLAB: How to divide the signal to AC and DC component.

After you calculate

Y = fft(signal,NFFT)/L;

then Y(1) is the average DC component, and Y(2:end) have to do with the AC components.

Best Answer

Related Solutions

MATLAB: I get the peaks at the wrong frequency and with the wrong amplitude.

MATLAB: How to divide the signal to AC and DC component.

Related Question