Need Help about FFT and STFT

M
Mikaela_sf · Jan 19, 2021 · 2.1K views
Question
I got a large amount of data, more than 5mb, so I can’t upload it. I would like to ask, my data is three xyz three-axis data collected from acceleratoter and it based on time, they are written in a .mat file. I want to intercept a small part of the time from this .mat file. What should I do? After that when i get the new data from that time ,how can i calculate the FFT Power Spectrum for a sliding window ? Can someone help me with this?
Expert Answer
Profile picture of Kshitij Singh
Kshitij Singh PhD Expert
Answered Nov 20, 2025
see the code below for (averaged) fft analysis and spectrogram
you can easily incorporate a test to select samples that matches the condition time > t_min and time < t_max.
%%%%

% load signal

%%%%


load('signal15.mat'% time (tout) and data (simout)

% t = tout;

channel = 1;

signal = simout(:,channel);

samples = length(signal);

Fs = 1/mean(diff(tout));


%%%%

% FFT parameters

%%%%

NFFT = 4096; %

OVERLAP = 0.75;

% 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;

%% decimate (if needed)

% NB : decim = 1 will do nothing (output = input)

decim = 2;

if decim>1

signal = decimate(signal,decim);

Fs = Fs/decim;

end

%%%%%

% 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
Have a different question? Ask here

Get a Free Consultation or a Sample Assignment Review!