Cannot perform FFT on a signal imported from Simulink

Base of my project is ADC set up in Simulink. I'm importing output volate of ADC from Simulink to Matlab workspace using "simout" block and "timeseries" format. I am able to plot this signal in time domain by the following:

plot(squeeze(V_out.data));

However I'm not able to get spectrum analysis of this signal using Fast-Fourier Transform (FFT) in Matlab. The following:

plot(fft(squeeze(V_out.data)));

results in "NaN + NaN * i".

If I try using "FFT" block in Simulink and then export this data to Simulink (via "simout" block), I get same values as for "V_out.data" except imaginary part is added, which is zeroed for all instances. Plotting this gives me flat function.

P.S.: Output signal of ADC is quanitized (sampled), if this matters for some reason. I'm also sharing ".mat" file of "V_out" ADC output voltage.

[Mathieu NOE]

hello

can you save your data (V_out.data) in a mat file and share it ?

[Luka Gacnik]

@Mathieu NOE I just edited the question with requested file.

Answer 1

[Mathieu NOE]

hello again

the issue was that the signal start with a NaN value

removing this first sample would solve all your issue. maybe you have to do something at your simulink model to avoid this NaN being created

do not pay attention to my bandpass filter code , this was already in the code ....

but I am surprised you set the sampling rate at 50 kHz for 1 Hz sine recording - quite overkill !

clc
clearvars

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = ('V_out.mat'); 
load(filename);
time = V_out.Time;% time vector 
signal = V_out.Data(:);
dt = mean(diff(time)); 
Fs = 1/dt;

[samples,channels] = size(signal);

% remove NaN values
indnan = find(isnan(signal));
signal(indnan) = [];
time(indnan) = [];

%% band pass  filter section %%%%%%
    f_low = 0.1;
    f_high = 100;
    N_bpf = 2;           % filter order
    [b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
    signal_filtered = filtfilt(b,a,signal);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000;    % 
OVERLAP = 0.75;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
subplot(211),plot(time,signal,'b');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / raw data ']);
xlabel('Time (s)');ylabel('Amplitude');
subplot(212),plot(time,signal_filtered,'r');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / filtered data ']);
xlabel('Time (s)');ylabel('Amplitude');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);

figure(2),plot(freq,20*log10(spectrum_raw),'b',freq,20*log10(spectrum_filtered),'r');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend('raw','filtered');

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 percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    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

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        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

[Luka Gacnik]

Thank you for such well explained answer. I don't know if you have any knowledge with ADCs or signal processing but shouldn't peak of frequency spectrum be near original signal's frequency? Because now it peaks towards DC (0 Hz). This is how V_out looks like (1 Hz input signal):

[Mathieu NOE]

hello

yes , I have some knowledge in these fields

simply here the sampling rate is super high (50 kHz) for a 1 hz sine wave

to see it in the fft plot, I had to decimate the data by factor 1000 to make it visible in the low frequency range

this is how the data and fft looks like after decimation by factor 1000

I also detrended the data before doing the fft

clc
clearvars

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = ('V_out.mat'); 
load(filename);
time = V_out.Time;% time vector 
signal = V_out.Data(:);
dt = mean(diff(time)); 
Fs = 1/dt;

[samples,channels] = size(signal);

% remove NaN values
indnan = find(isnan(signal));
signal(indnan) = [];
time(indnan) = [];

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1000;
if decim>1
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    Fs = Fs/decim;
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 250;    % 
OVERLAP = 0.75;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
plot(time,signal,'b');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / raw data ']);
xlabel('Time (s)');ylabel('Amplitude');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(detrend(signal),Fs,NFFT,OVERLAP);

figure(2),plot(freq,20*log10(spectrum_raw),'b');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');

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 percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    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*ones(1,channels));
        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

[Luka Gacnik]

@Mathieu NOE I apologize in advance for asking additional question, but what is decimation about? What does it do and why? Otherwise again well explained, I couldn't have asked for more :)

[Mathieu NOE]

hello Luka

no problem - normally we would define the sampling rate as at least 2 x the highest analog signal frequency to be sampled (Shannon thorem). in practice we may take a larger factor like 5 or 10;

of course you could do a much faster sampling but what for ?

you have to remember the relationships between signal sampling rate , duration and resulting fft proprerties :

• the fft frequency vector goes from 0 to Fs/2 (Nyquist frequency)
• the fft frequency resolution (delta f , or df in short) is related to Fs and the buffer length (fft length = nfft) : df = Fs/nfft.
• so signal duration, sampling rate will impact what we achieve on the fft results

in your case, in order to distinguish the 1 Hz peak you need a delta f (df) lower than 1 (says 0.5 or 0.25 Hz => 1/df = 2 or 4 here) ;

if you have a high Fs , this means you need to record a very long signal of length > 2 or 4 times Fs (this is large !!) , and also you do a huge fft computation on this large data to finally only use the first fft points (bins) and throw 99% to the bin.

that is alot of work for nothing

better downsample (in matlab words : decimate) so that Fs is much lower , which is anyway not an issue as we are looking for low frequency fft data.

this shorten the data length to be fft processed and save a lot of computation burden

see also FYI :

Decimation - dspGuru

Downsampling (signal processing) - Wikipedia