2
$\begingroup$

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
$\endgroup$
2

2 Answers 2

2
$\begingroup$
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.

$\endgroup$
4
  • $\begingroup$ That feels a bit simplistic: the correct normalization depends on the nature of the signals and for broadband signals its the RMS not the mean. $\endgroup$
    – Hilmar
    Aug 14, 2023 at 12:22
  • $\begingroup$ edited with a link to Dan's answer on the same subject $\endgroup$
    – Jdip
    Aug 14, 2023 at 12:27
  • $\begingroup$ Adding answer from Jdip to my fft function gives me the expected results, thanks a lot. Follow-up question based on Hilmar and Dan Boschen: Does that mean as long as I perform ffts only for signals with distinct frequency components I am good to go using the mean, but if I wanted to perform an fft of for example a white noise signal to check if it has the expected bandwidth I would have to use RMS? $\endgroup$ Aug 14, 2023 at 12:38
  • 1
    $\begingroup$ That’s correct Gerald. See this answer as well. Don’t forget to mark the answer to close it. $\endgroup$
    – Jdip
    Aug 14, 2023 at 12:48
2
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Certainly, this is not a comment; it is a full-fledged answer. In addition to being to the point and useful in a wider application area, the answer offers a relevant code review. $\endgroup$
    – V.V.T
    Aug 15, 2023 at 5:24
  • $\begingroup$ I agree with V.V.T. Hadn't heard of window overlap before. My signal's frequency spectrum is not changing over time, but do you think using overlap would improve my measurement quality anyway just because I can average the fft results over more windows? $\endgroup$ Aug 15, 2023 at 6:22
  • 1
    $\begingroup$ Hi Geralt, in general overlapping windows are useful because they allow more fine-grained views of signal changes over time, or a more "correct" average result of the whole signal, by averaging results across windows, if that is desired instead. I think the answer is yes, especially if there is noise in your actual data. $\endgroup$
    – user25849
    Aug 16, 2023 at 3:19

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.