MATLAB Answers

Need Help about FFT and STFT

58 views (last 30 days)
Xiaoming Hu
Xiaoming Hu on 12 Jan 2021
Edited: Xiaoming Hu on 22 Jan 2021
Hi!
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? Best Regards!

  0 Comments

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 13 Jan 2021
hello
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

  5 Comments

Show 2 older comments
Xiaoming Hu
Xiaoming Hu on 15 Jan 2021
thanks ! one more question:
will the window lengths and overlaps influences the outcome ? How can i select a good size of window?
Mathieu NOE
Mathieu NOE on 18 Jan 2021
hello
If your window lenghth is nfft , the spectrum frequency resolution is df = Fs/nfft ; so there are compromise between frequency resolution and time resolution for a spectrogram plot.
I usually start with an overlap of 50 to 75% of nfft .
if I really need to heave a smoother and refined time resolution , I can go up to 95 / 97% of overlap.
Xiaoming Hu
Xiaoming Hu on 22 Jan 2021
Thanks Professor!!! HHH,
I have finisch the nfft and i just get the area what i want to use from 0.9hz to 15hz. Then i want to get a root mean square value. The code what i have did till right now is below.
windowsize = Fs*2; %because the window is 2s , and sample frequncy is 4096, so it is the size
window = hamming(windowsize);
nfft = windowsize; %equal to window size
noverlap = Fs; % overlap is 1s,so it is Fs
xR2 = ACC_DATA_R(3,time_postural_tremor);%tremor data is z axis used
[S2,F2,T2,P2] = spectrogram(xR2,window,noverlap,nfft,Fs);%%this is for calculate the average power
f_bigger_09Hz = find(F2>0.9,1,'first');
f_15Hz = find(F2<=15,1,'last');
%Right
figure;
plot(F2(f_bigger_09Hz:f_15Hz),P2(f_bigger_09Hz:f_15Hz,:))
ylabel('power')
xlabel('Frequency in Hz right')
grid on
%Left
figure;
plot(F1(f_bigger_09Hz:f_15Hz),P1(f_bigger_09Hz:f_15Hz,:))
ylabel('power')
xlabel('Frequency in Hz left')
grid on %% iget the area what i need(0.9hz to 15h
Can u give me a help how to calculate this area to get a root mean square value?
Thanks!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by