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