2
$\begingroup$

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')
$\endgroup$
8
  • $\begingroup$ This is in my opinion the best reference on the subject. $\endgroup$
    – Jdip
    Nov 26, 2022 at 13:17
  • 1
    $\begingroup$ Yes, it is. But I don't find an answer to my question in it. Thank you again btw $\endgroup$
    – user49811
    Nov 26, 2022 at 16:25
  • 1
    $\begingroup$ Ok... In Matlab documentation is named "Hanning". I didnt' know that. Do you have any idea about my problem $\endgroup$
    – user49811
    Nov 29, 2022 at 14:10
  • 1
    $\begingroup$ Can you share your code? So we can have something to work with ;) $\endgroup$
    – Jdip
    Nov 29, 2022 at 23:25
  • 1
    $\begingroup$ @Jdip I attached a sample one $\endgroup$
    – user49811
    Nov 30, 2022 at 6:51

1 Answer 1

1
$\begingroup$

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.

enter image description here

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. enter image description here

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$: enter image description here

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])
$\endgroup$
2
  • $\begingroup$ I am not interested in the PSD. I just want a spectrum with definition similar to the average fft (without all that peaks that make it difficult to read), but without an upward shift. Is that impossible? $\endgroup$
    – user49811
    Nov 30, 2022 at 12:57
  • $\begingroup$ @user49811 let's take this to a chat: chat.stackexchange.com/rooms/140951/fft-shift $\endgroup$
    – Jdip
    Nov 30, 2022 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.