Overlapping windows using the fft for time discrete data

Hello,

I want to calculate the frequency domain of a measured signal by using the fft. My measurements have the samplerate of 5000Hz over 45s. The signals are acceleration in x,y,z direction of my coordinate system. I need the frequency domain signal for every second. The window I use is the Hann-Window and I decided to use the overlap of 50% as usual.

The problem is:

• if I say the win_length is 5000Hz as the samplerate, then the result is for every half second
• if I say the win_length is 10.000Hz then the result is for every second, but the algorithm cosiders 2s.

I must be wrong in understand the fft for a time discrete signal, but I read a lot about that. So what am I doing wrong?

sens_length = length(sensordata);
w = hann(win_length);
num_freq = win_length/2 ;
num_w = 2*sens_length/win_length-1;


for t=1:num_axis % for every axis of the coordinate system
    for i = 1:num_w   %for every window
        sensor_freq(:,i,t) = fft(w .* sensordata((i-1)*num_freq+1:(i-1)*num_freq+win_length,t));
        sensor_freq_abs(:,i,t) = abs(sensor_freq(1:num_freq,i,t));
    end
end

This is how one of my signals after measurement looks like:

[Mathieu NOE]

Hi

what you are doing in the inner loop do already exist (in better coding ) with matlab's pwelch function

so you can reduce the complexity (and remove the bug)

also I se your signal is quite unsationnary

consider doing time / frequency analysis using specgram / spectrogram

[AC_GBX]

Hi,

what do you mean by "time / frequency analysis using specgram / spectrogram"?

I tried pwelch like this:

[Y,f] = pwelch(sig_act(2:255001,1,1),hann(5000),2500,5000,5000);

The problem is: For all the measured 45s I only get one vector if I use pwelch like this. But I want to get a vector with the amplitudes for every second. So in total 45 vectors containing the amplitudes for the frequency from 1 to 2500Hz.

[Mathieu NOE]

so what is the purpose of the 45 spectrum vectors ? see evolution of spectrum with time ?

this is exactlty what spectrograms are intended for. There are already functions and demos

see the signal processing doc : Practical Introduction to Time Frequency Analysis

[AC_GBX]

Hi Mathieu,

thank you for your answer. The purpose of the 45 spectrum vectors is that I look at 45 operating states and I assume that one second represents one operating state.

The function spectrogram nearly does what I want. But I am not sure if I use it the correct way, because my results are not what needed.

I used it like this:

f = linspace(1,5000,5000)';
[x,w,t] = spectrogram(GBX_LH_act_x(1:225001),hann(5000),2500,f,5000);
x_oneside = abs(x(1:2500,:)); %only the amplitude for one side

But the result is very similiar to my previous calculation only using the fft.

I get in total 89 vectors if I use a windowlength of 5000 and I get a 44 vectors for windowlength of 10000.

What am I doing wrong?

I was asking myself I it would be better to overlap the fft during one second and calculate the mean of the windowed results just like pwelch does it. Do I need to do that?

[Mathieu NOE]

I prefer to "redesign" your code according to my own ideas... hope it matches your needs

so I ended up doing 2 things :

• one example how to use pwelch on a buffer of 4096 samples, on each I do a pwelch averaging fft (with nfft = 512). That gives me 62 spectra that I display on a graph
• second example is based on spectrogram. It's very similar in naturen bt data are displayed as time / frequency / amplitude (color coded , levels are in dB). It uses the funtion : spectrog_peak.m (attached)

hope this helps you .

% dummy signal
Fs = 5000;
samples = 255001;
dt = 1/Fs;

t = (0:dt:(samples-1)*dt);
omega = 2*pi*(25+20*t);

sensordata = randn(size(t)) + 5*sin(omega.*t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

samples = length(sensordata);

num_axis = 1;

NFFT = 512;    % to have highest frequency resolution , NFFT should be greater than typical max length of files ( 1 second duration at Fs = 16 kHz = 1600 samples);
NOVERLAP = round(0.75*NFFT);

samples_per_frame = NFFT*8;

num_w = floor(samples/samples_per_frame);

% your code modified - multiple spectra 

w = hamming(NFFT);

sensor_spectrum_all = [];

for t=1:num_axis % for every axis of the coordinate system
    for ci = 1:num_w   %for every window
        ind_start = 1+(ci-1)*samples_per_frame;
        ind_stop = ind_start + samples_per_frame-1;
        [sensor_spectrum, freq] = pwelch(sensordata(ind_start:ind_stop),w,NOVERLAP,NFFT,Fs);
        % concatenate spectra
        sensor_spectrum_all = [sensor_spectrum_all sensor_spectrum];
        
    end
end

figure(1),plot(freq,20*log10(sensor_spectrum_all));
title(['Spectra  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz  / Data buffer ' num2str(samples_per_frame) ' samples' ]);
xlabel('Time (s)');ylabel('Frequency (Hz)');

% spectrogram demo
Satur = 80;  % dB range scale (means , the lowes displayed level is 80 dB below the max level)
fmin = 1;
fmax = Fs/2;

% option_fw = 1; % FFT pond L
% option_fw = 2; % FFT pond A
option_fw = 1;
tmin = min(t);    tmax = max(t);
[sg,fsg,tsg] = spectrog_peak(sensordata,NFFT,Fs,0.75,Satur,tmin,tmax,fmin,fmax,option_fw);

% plots sg
figure(2);
imagesc(tsg,fsg,sg);axis('xy');colorbar('vert');
title(['Spectrogram  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

[AC_GBX]

Thank you very much. This is very helpful and might work for my problem.

Is there any possibility that I receive a vector with a result for every frequence? The vector freq must be 5000x1. My result sensor_spectrum must be the same length as the observed frequence (5000Hz). I thought normally this can be made with the variable f in pwelch, but if nfft is used there seem to be no possibility to use f as well.

[Mathieu NOE]

hello again

some remarks :

1. when you do an fft (as pwelch do) there is a relationship between the frequency resolution (df) with Fs (500 Hz) and NFFT : df = Fs/NFFT so parameters are linked.
2. the choice of NFFT (and overlap) is a compromise (your choice) between frequency resolution (higher NFFT) and lot of averaging (smaller NFFT and maximum OVERLAP) if the data are corrupted with noise. You change change the NFFT and samples_per_frame values to see the effect on your data
3. if Fs = 5000 Hz, the maximum observable frequency you can plot is Fs/2 = 2500 Hz. This is Shannon theorem.
4. if you have a need to have the ouput spectrum length equal to 5000 , you can always interpolate on a new frequency vector , example below :

freq_interp = linspace(1,Fs/2,5000); % create a new freq vector of length 5000 
sensor_spectrum_interp = interp1(freq,sensor_spectrum,freq_interp)

you can insert these lines before doing the spectra concatenation

does it answer your question ?

[AC_GBX]

Yes this answers my questions. Thank you very much. You helped a lot. I will try to interpolate the first half of the vector up to 2500 values due to shannon .

[Mathieu NOE]

ok

FYI, the freq vector generated from the pwelch function is already in range 0 to Fs/2

maybe there is no need to take again the first half of the freq vector (or you will be then limited to 0 : Fs/4 range )