Amplitude scaling of window functions for FFT

Question

I need to perform an FFT on a signal sampled with 20 kHz and a measurement time of about 10 seconds. The signal contains frequencies of up to 2 kHz but I am mainly interested in the bandwidth of 0 to 300 Hz.

Just applying Matlab's fft function gives me an unnecessarily high frequency resolution that results in a rather noisy plot so I started looking into windowing. The idea is to split the unnecessarily long measurement into multiple windows of defined length, applying the fft function to each of these windows, and afterwards average over all windows.

To avoid leakage effects, I applied a Blackman-Harris window function to each of the windows before applying the fft function. According to National Instruments documentation, it is a good general-purpose function with good amplitude accuracy due to its wider main lobe compared to uniform or Hann window.

Before applying my script to the real data, I thought of generating some artificial data to test the behavior of my script and choose the right window size and function. My first observation was that the Blackman-Harris window requires a much larger window size than the uniform window to be able to tell apart two close frequencies. This is okay and explainable for me, as the Blackman-Harris main lobe is relatively wide. My second observation however was that the amplitudes all seem to be off by roughly the same factor compared to the fft with the uniform window.

So my questions are: Does this amplitude scaling result from the window function scaling down the measurement data points in the time domain? And can this scaling factor be reliably calculated and added to the fft?

Below you can find the Matlab script that I used so you can reproduce my findings. Signal y1 is just for reference. Signals y2 and y3 should test how good close frequencies can be distinguished. Signal y4 should test how well the window function handles leakage effects.

x = 0:5e-05:10; % 10 seconds of measurement at 20 kHz
y1 = sin(x * 2 * pi * 30.0 + 0.9) * 1.0; % 30.0 Hz, AMP 1.0, Phase 0.9
y2 = sin(x * 2 * pi * 75.0 + 1.5) * 0.8; % 75.0 Hz, AMP 0.8, Phase 1.5
y3 = sin(x * 2 * pi * 77.0 + 2.6) * 0.2; % 77.0 Hz, AMP 0.2, Phase 2.6
y4 = sin(x * 2 * pi * 91.3 + 3.1) * 0.5; % 91.3 Hz, AMP 0.5, Phase 3.1
y = y1 + y2 + y3 + y4;

% Variation of window size
ws = [5000, 20000, 40000, 120000];

figure;
for ii = 1:numel(ws)
    [f, P] = myfft(y, 20000, ws(ii));
    plot(f, P, 'DisplayName', ['Signal length: ' num2str(numel(x)) '; Window size: ' num2str(ws(ii))]); hold on;
end
title('Variation of window size');
xlabel('Frequency in Hz');
ylabel('Absolute amplitude');
xlim([0 120]);
grid on;
legend;

This is my custom FFT function. Under the "Perform transformation" line you can change the w to a 1 in order to get results with uniform window.

function [f, P] = myfft(data, fs, windowSize)

    % Determine singal length
    L = length(data);

    % Ensure data is a column vector
    if size(data,1) < size(data,2)
        data = data';
    end

    % Create weighting function
    w = blackmanharris(windowSize);

    % Init result averaging vector
    P_ii = [];

    % Iterate through windows
    for ii = 1:windowSize:L
        
        % Only process full-length windows
        if L - ii >= windowSize - 1

            % Perform transformation
            Y = fft(data(ii:(ii+windowSize-1)) .* w);
            P2 = abs(Y/windowSize);
            P1 = P2(1:windowSize/2+1);
            P1(2:end-1) = 2 * P1(2:end-1);

            % Deal with row / column vector issues
            if size(P1,1) > size(P1,2)
                P_ii = [P_ii, P1];
            else
                P_ii = [P_ii, P1'];
            end

        end

    end

    % Init result vector
    P = zeros(windowSize/2+1,1);

    % Averaging
    for ii = 1:size(P_ii,1)
        P(ii,1) = mean(P_ii(ii,:));
    end

    % Frequency band
    f = fs*(0:(windowSize/2))/windowSize;
    f = f';

end

Answer 1

P2 = abs(Y/windowSize);

Here you are scaling your DFT coefficients by the window length. This is correct for a uniform (rectangular) window. However, to have a consistent value for the DFT coefficients regardless of the window weights $w[n]$, the correct normalization factor is $$S = \sum_{n=0}^{N-1}w_n$$ This normalization coefficient takes into account both the length of the window $N$ and the window gains $w[n]$ (note that for a uniform, i.e. rectangular window, $S=N$).

You should divide each coefficient by this factor to retain energy consistency across different windows of different lengths.

For more details (addressing Hilmar's comment), see this answer by Dan Boschen and this one as well.

Answer 2

I don't have enough rep for comment but simply wanted to point out that your code yields windows with zero overlap, which would leave big gaps in your analysis were your signal to change with time (perhaps it doesn't). For your reference, the Blackman-Harris window achieves constant-overlap-add at 2/3 overlap | 1/3 hop size. So perhaps instead of

for ii = 1:windowSize:L

doing

for ii = 1:round(windowSize/3):L 

would be good. Also, if you're looking for something with less leakage than a uniform window but higher frequency resolution, a simple raised cosine (Hann) window is always good, and does constant-overlap-add at a much more convenient 1/2 | 1/4 | 1/8 etc. hop size.

I looked up the relevant question about comment vs answer in meta.stackexchange and it appears that people with <50 rep are actually expected to post comments as answers. Sorry if I'm wrong.