Cannot perform FFT on a signal imported from Simulink

4 views (last 30 days)
Luka Gacnik
Luka Gacnik on 11 Apr 2022
Commented: Mathieu NOE on 11 Apr 2022
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.

Answers (1)

Mathieu NOE
Mathieu NOE on 11 Apr 2022
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
  4 Comments
Mathieu NOE
Mathieu NOE on 11 Apr 2022
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 :

Sign in to comment.

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!