upward shift in FFT windowed - attenuate fluctuation in FFT

Question

I would like to limit the noise after the FFT. I have used this strategy:

1. Divide the signal in segments using a Hann window
2. Perform abs(FFT(segment))
3. Ensemble average
4. Apply the gain factor of the Hann window.

This ends up being the same as doing the square root of the PSD obtained with the Welch method. I am satisfied in terms of reduced fluctuation with this "filtered" spectrum, but it seems to be shifted above, and in my case the magnitude of the spectra is important.

I found a similar behavior of what I am saying in Fig. 2 here.
Is there a way of avoiding or compensating this upward shift? I have noticed that avoiding the abs at step 2 results in avoiding this upward shift, but it causes also more fluctuations, which makes the filtration useless. Are there other method to reduce fluctuations? Here it is a sample code

clear; clc; close all;

rng default

n = 1:10000;

Sn=0.1;
Fs=1/Sn;
L=length(n)

duration = 1000;

x = pinknoise(duration*Fs);

freq_original = transpose(Fs*(0:(L)/2)/L);
W_len=500;
freq = transpose(Fs*(0:(W_len)/2)/W_len);

%Original Spectrum
xf_original=abs(fft(x))./L
xf_original= xf_original(1:L/2+1,:);
xf_original(2:end-1,:) = 2*xf_original(2:end-1,:);


%Division of the signal in segment and Henning windowing
xdiv = buffer(x(:,:),W_len,(W_len)/2, 'nodelay')
xdiv=xdiv(:,1:38)
A = hanning(W_len); 
xdivHanning=A.*xdiv


%fft  and absolute value of every window
xf2=abs(fft(xdivHanning))./W_len
xf2= xf2(1:W_len/2+1,:);
xf2(2:end-1,:) = 2*xf2(2:end-1,:);
xf2=transpose(mean(transpose(xf2)))


%Amplitude Gain Factor for Hanning window
CrrFac=2

figure('DefaultAxesFontSize',13)
loglog(freq_original ,(xf_original),'r')
hold on
plot(freq ,CrrFac.*abs(xf2),'b')

legend('original spectrum','Upward shifted fft')

Answer 1

The paper you linked to explains why there is a shift. It has to do with the length of the input data, which is directly related to frequency resolution and the distribution of noise power in each fft bin. With larger windows, frequency resolution $d_f$ gets smaller (better!) and noise gets spread out.

It all comes down to the question: what are you interested in? I'm assuming your signal is somewhat stationnary (otherwise you'd need time-frequency analysis)

• If you're interested in exact signal amplitudes for a discrete set of frequencies, compute a straight fft (choose the length according to what frequencies you're interested in) and read the amplitudes out.
  If you're not sure which frequencies you're interested in, for example you want to find the 3 dominant frequencies of a signal, you can use frequency averaging as long as you have sufficient frequency resolution (choose the window length accordingly). In both cases, you'll need to scale according to the window used.
• If you're interested in noise power, you should use the Power Spectrum or Power Spectral Density and use the appropriate scaling factors as well.

---

For example, I added two sine waves (amplitude $A = 3$, frequencies $0.5\,\text{Hz}$ and $0.6\,\text{Hz}$) to pink noise, with $F_s = 10\,\text{Hz}$ and $N = 10000$. I'll compare doing a straight fft ($N=10000$ so $d_f = 0.001\,\text{Hz}$) vs averaging a bunch of length-500 ffts together ($d_f = 0.02\,\text{Hz}$).

Here are the results if you're interested in the signal amplitudes. Notice the amplitudes are correct, but the noise floor is "shifted" when frequency averaging.

Here are the results if you're interested in noise. Notice the noise power is consistent with each method, but the signal power magnitudes are now shifted when doing frequency averaging.

Finally, here is where frequency resolution matters: change the sine waves to $0.5\,\text{Hz}$ and $0.51 \,\text{Hz}$, the difference is now below $d_f$:

Here is Matlab code to generate these:

    clear; clc; close all;
rng default

% test signal
n = 1:10000;
Sn = 0.1;
Fs = 1/Sn;
L = length(n);
duration = 1000;
x = pinknoise(duration*Fs);
sig = 3*(sin(2*pi*0.5/Fs*n)' + sin(2*pi*0./Fs*n)');
x = x+sig;


% Straight fft and PSD with magnitude squared method
xf = abs(fft(x))/L;
xf= 2*xf(1:L/2+1,:);

xf_original = abs(fft(x)).^2;
xf_original= xf_original(1:L/2+1,:);
xf_psd = 2*xf_original/(Fs*L); % omit Fs for Power Spectrum

% frequency averaging, and PSD with welch's method
W_len = 500;
xdiv = buffer(x(:,:), W_len, W_len/2, 'nodelay');
A = hanning(W_len); 
xdivHanning = A .* xdiv;

xf2_avg = mean(abs(fft(xdivHanning)/sum(A)),2);
xf2_avg = 2*xf2_avg(1:W_len/2+1,:);

xf2 = abs(fft(xdivHanning)).^2;
xf2 = xf2(1:W_len/2+1,:);
xf2 = transpose(mean(transpose(xf2)));
xf2_psd = 2*xf2/(Fs*sum(A.^2)); % omit Fs and change to sum(A).^2 for Power Spectrum

p = pwelch(x, W_len, W_len/2, W_len, Fs);

% frequency vectors
freq_original = transpose(Fs*(0:(L)/2)/L);
freq = transpose(Fs*(0:(W_len)/2)/W_len);

% Plot
figure(1)
plot(freq_original, 10*log10(xf_psd))
hold on
plot(freq , 10*log10(xf2_psd))
plot(freq , 10*log10(p))
ylim([-30 50])
legend('original PSD', 'Welch PSD', 'Matlab Welch PSD')
title('PSD')
grid on

figure(2)
subplot 211
plot(freq_original, xf)
hold on
plot(freq , xf2_avg)
grid on
title('Get frequency amplitudes')
legend('straight fft', 'average fft')
subplot 212
plot(freq_original, xf)
hold on
plot(freq , xf2_avg)
grid on
xlim([0.3 0.8])