Differences between FFT methods

Hello everyone, I have been trying some different methods for computing the FFT of some accelerometer data I collected. The first method is a regular FFT with an Hamming window across all samples. The second and third methods are the STFT and welch method with the parameters that you can see below. Since these methods have differences in the way they perform the calculations, I was expecting some variation of the results. However, it appears that these are two big. I have to present these results and want to know which is the best method. I was wondering if anyone could help me better understand what is actually happening that causes these differences:

1 - Should I be diving the spectra by sum(w), honestly just saw that it was a good practice in FFT calculations so I am doing it. But to be honest it doesn't make much sense to me to present the results in -dB, so I'm not so sure about that. SInce this is vibration data wouldn't it make more sense that the peaks are, for example, 60dB, 100dB and so on? If I don't normalize the results (divide by sum(w)) the amplitudes increase but the pwelch method still has values below 0dB, which I don't understand as well.

2 - My signal is stationary, does it make sense to be using the STFT in such situation since I don't need time information?

3 - How could I average the results from the STFT or FFT to get a more smooth line like in the pwelch method?

I am trying to fully understand what happens in the background of these calculations and would be glad to have some insights from more experienced people on the topic. Thank you very much in advance.

% Signal Parameters
fs_m = 26700;            % Sampling Frequency (MEMS)
Nsamples1 = 775862;      % Number of Samples

% Window parameters 
wlen = 2^14;                                   % Window length (recommended to be power of two)
nfft = 2*wlen;                                 % Number of fft points 
overlap = 0.50;                                % Defining overlap as 75%
hop = wlen * (1-overlap);                      % Hop size 
w = hanning(wlen, 'periodic');                 % Hanning window

% Testing pwelch method [Pxx,F] = pwelch(X,WINDOW,NOVERLAP,NFFT,Fs) 

[Sxx,f_welch] = pwelch(acc_MEMS1, w, wlen-hop, nfft, fs_m); % single-sided PSD
abs_Sxx = abs(Sxx);
abs_mag = 2*abs_Sxx(1:nfft/2+1);
SxxdB = 20*log10(abs_Sxx);

% Perform STFT
[s, f, t] = stft(acc_MEMS1, fs_m, "Window", w, "OverlapLength", wlen-hop, "FFTLength", nfft);

% Filter out frequencies less than 0
positive_f = f >= 0; f_pos = f(positive_f,:); s_positive = s(positive_f, :);
sf_w = s_positive/sum(w);
% sf_w = s_positive;
abs_sf = abs(sf_w); avg_mag_s = 2*abs_sf(1:nfft/2+1); avg_mag_sdB = 20*log10(avg_mag_s);

% FFT
[fft_y_Hann, fft_y_AccH] = fft_analysis_Hann(Nsamples1, timeVector1, acc_MEMS1, fs_m);
%% Inside the function fft_analysis_Hann

if mod(N,2)~=0 % da Berechnungsformeln nur mit N/2 integer ohne Warnung funktioniren (letzten Wert abschneiden)
    N = N-1;
    t = t(1:end-1);
    y = y(1:end-1);
end

w = hanning(N);  % Hanning window
f_accH =(0:(N/2))*f_s/N; % only frequency >0

y_window = y(:).*w(:);
fft_y_komp = fft(y_window)/sum(w); % Normalization of the Amplitude
% fft_y_komp = fft(y_window); % No normalization

abs_hwindow = abs(fft_y_komp);     %  Window
fft_y_Hann = 2*abs_hwindow(1:N/2+1);

%% End of the function

fft_y_HanndB = 20*log10(fft_y_Hann);

% FREQUENCIES OF INTEREST
fdesM = 3000; fdesP = fdesM;

%Plot

f1 = figure('Name', 'Measurement Results');   

plot(fft_y_AccH, fft_y_HanndB, 'Color', [0.169 0.169 0.169 0.8], 'LineStyle', ':', 'LineWidth', 0.8); xlabel("Frequency"); ylabel("Amplitude (dB)"); pbaspect([4 1 1])
hold on;
plot(f_pos, avg_mag_sdB, 'Color', [1 0 0 0.6], 'LineStyle', '-', 'LineWidth', 1.2); xlabel("Frequency (Hz"); ylabel("Amplitude (dB)"); title("Frequency Domain Comparison Measurement N°1"); xticks(0:60:fdesM); xlim([0 fdesP]);  
pbaspect([4 1 1]);
plot(f_welch, SxxdB, 'Color', [0 0 1 0.6], 'LineStyle', '-', 'LineWidth', 1.2); xlabel("Frequency (Hz"); ylabel("Amplitude (dB)"); title("Frequency Domain Comparison Measurement N°1"); xticks(0:60:fdesM); xlim([0 fdesP]);
pbaspect([4 1 1]); legend('Regular FFT', 'STFT', 'Welch Method')
hold off;

Answer 1

[Sahaj]

Hi Ivo. Here are the answers to your questions:

1. Dividing the spectra by the sum of the window function (sum(w)) is a common practice in FFT calculations. This normalization is done to ensure that the amplitude values are consistent across different window lengths and types. By normalizing, the amplitudes are scaled relative to the window function, which helps in comparing different spectra. However, the choice of presenting the results in dB is subjective and depends on the context. If you prefer to present the amplitudes in linear units, you can skip the normalization step. The negative dB values you observe in the pwelch method are a result of the logarithmic scaling used in the conversion to dB. In the logarithmic scale, values below 0 dB are represented as negative dB values. It's important to note that the negative dB values do not indicate negative amplitudes but rather lower amplitudes relative to the reference level.
2. The Short-Time Fourier Transform (STFT) is commonly used for analyzing non-stationary signals where time information is important. Since you mentioned that your signal is stationary and you don't require time information, using STFT may not be necessary in your case. The STFT provides a time-frequency representation by dividing the signal into short segments and calculating the FFT for each segment. If you only need the frequency domain information, using a regular FFT may be more appropriate and computationally efficient.
3. To obtain a smoother line in the frequency domain, you can average the results from the STFT or FFT. This averaging can be done by taking multiple segments of the signal and calculating the FFT or STFT for each segment. Then, you can average the spectra across these segments to obtain a smoother representation. Averaging helps to reduce the effects of random variations and noise, providing a more reliable estimate of the underlying frequency content.

Regarding what is happenning in the background of these calculations:

Pwelch Method:

1. The pwelch function is called with the input signal acc_MEMS1, the Hanning window w, the overlap length (wlen-hop), the number of FFT points nfft, and the sampling frequency fs_m.
2. The output includes the power spectral density (Sxx) and the corresponding frequency values (f_welch).
3. The absolute value of Sxx is calculated (abs_Sxx).
4. The magnitudes of the positive frequencies up to nfft/2+1 are extracted (abs_mag).
5. The magnitudes are converted to dB using the formula 20*log10(abs_Sxx) to obtain SxxdB.

STFT Method:

1. The stft function is called with the input signal acc_MEMS1, the sampling frequency fs_m, and additional parameters including the Hanning window w, the overlap length (wlen-hop), and the number of FFT points nfft.
2. The STFT function divides the input signal into overlapping segments of length wlen with the specified overlap.
3. The Hanning window is applied to each segment.
4. The FFT is performed on each windowed segment to obtain the frequency spectrum for that segment.
5. The resulting spectra for each segment are stored in s.
6. The corresponding frequency values are stored in f, and the time values are stored in t.

Regular FFT with Hanning Window:

1. The function fft_analysis_Hann is called with parameters including the number of samples (Nsamples1), the time vector (timeVector1), the input signal acc_MEMS1, and the sampling frequency fs_m.
2. Inside the function, the number of samples is checked to ensure it is even. If it is not, the last sample is removed.
3. A Hanning window of length N is created.
4. The frequency values for the FFT output are calculated using f_accH.
5. The input signal is multiplied by the Hanning window to apply the windowing.
6. The FFT is performed on the windowed signal, and the resulting spectrum is stored in fft_y_komp.
7. The amplitude spectrum is normalized by dividing it by the sum of the Hanning window.
8. The absolute value of fft_y_komp is calculated to obtain the windowed spectrum (abs_hwindow).
9. The magnitudes of the positive frequencies up to N/2+1 are extracted and stored in fft_y_Hann.
10. The magnitudes are converted to dB using the formula 20*log10(fft_y_Hann) to obtain fft_y_HanndB.

[Ivo Bastos]

Thank you for your answer! Suppose I would like to calculate the fft of a sine chirp that goes from 1-5000Hz. Since this signal's frequency changes over time, I should use either the STFT or the pwelch method correct? Also, averaging the STFT/FFT results is pretty much the same as just using the pwelch function?

Based on what you mention regarding the dB scale, in this case, with normalization, all the values are below 0dB. Why exactly are all the amplitudes lower than the reference level? That does not seem to make sense to me. I am sorry if it is a ridiculuous question but I trully don't understand