FFT of Large CSV File - Fundamental doesn't match data or simple sine test case

I'm trying to take the FFT of a very large .csv file from and oscilloscope, but the results do not correlate to what the oscilloscope says, nor a sample test case (included in the below script). I've checked my FFT script against several online examples and it seems to be correct. Any help is appreciated.

SampleInt = 0.00000004; % interval at which scope samples
ScopeFreq = 16.667; %[Hz] Reading from O-scope

% Spreadsheet Info:
spreadsheetname = 'tek0001ALL.csv';   %If it's in the current directory... need to conver to .xlsx file.

%Define Variables for Matlab
ScopeTime = csvread(spreadsheetname,475298,0,[475298,0,1974610,0]);
disp('Scope Time Imported, Script Line 23')
Va = csvread(spreadsheetname,475298,1,[475298,1,1974610,1]);
disp('Scope Va Imported, Script Line 25')
disp('Excel Importing Done!!!')

%Define Time Domain Series for Plotting
timelength = 0:1:(length(Va)-1);
Time = timelength.*SampleInt;

%Calc' the RMS value of Va
Va_RMS =rms(Va)
Max = max(Va)

% Do the FFT!!
L=length(Va)
Dt = SampleInt;        % Get the sampling rate
Fs = 1/Dt; 
NFFT = 2^nextpow2(L);  % automaticlaly figures out what n is for 2^n
X=fft(Va,NFFT/2);
mx = abs(X);
f =  (0:NFFT/2-1)*Fs/NFFT; 

%Plot just the FFT
figure(1)
stem(f,mx)
title('FFT of Raw Test Data')
xlabel('Freq. [Hz]')
xlim([0,200])

%Now plot Va against a sinewave to check fundamental from FFT
Vv=Max*sin(2*pi*16.8*Time);
figure(2)
plot(Time,Va,Time,Vv)
xlabel('Time [Sec]')
grid on

% Now Check the FFT function using the reference sinewave
Lv=length(Vv)
Dtv = SampleInt;        % Get the sampling rate
Fsv = 1/Dtv; 
NFFTv = 2^nextpow2(Lv);  % automaticlaly figures out what n is for 2^n
Xv=fft(Vv,NFFTv/2);
mxv = abs(Xv);
fv =  (0:NFFTv/2-1)*Fsv/NFFTv; 

figure(3)
stem(fv,mxv)
title('FFT of Raw Test Data')
xlabel('Freq. [Hz]')
xlim([0,200])

Note that when I run this script the FFT produced a fundamental frequency of ~ 12Hz, but that the fundamental is actually 16.67Hz.

Answer 1

[David Goodmanson]

Hi Calvin,

this is yet another example of the destructive power of using nextpow2 with the fft in the wrong circumstances.

I think the 2^n thing is a relic from 50 years ago. Under the right circumstances, say avoiding a large prime number of fft points, it makes some sense. Most of the time it's unneeded, and often it's counterproductive.

For the frequency grid, with the nextpow2 method the correct version of fv would be

fv =  (0:NFFTv/2-1)*Fsv/(NFFTv/2); 

That's because the spacing for an fft frequency array is Fs/npoints and you have npoints = NFFTv/2. However, your original signal is a very nice single oscillation of a sine wave in 1499312 points. If you fft that, you get a peak at the correct frequency. When you truncate the wave by using nextpow2/2 as you did, or pad it with zeros using nextpow2, you no longer have a nice continuous oscillation. Two things happen, both bad:

[1] After the continuous oscillation is messed with, the original single frequency peak picks up multiple frequency components, appropriate to describe the cutoff wave. But that is not what you want.

[2] the accompanying frequency grid no longer has the appropriate spacing to represent the original sine wave frequency. In your case you get a peak at the closest available frequency and some spurious confetti due to truncation at a bunch of other frequencies, including DC.

When npoints is a large prime number the fft is slow, and in that case it's not a bad idea to alter the number of points. In your case the fft is so fast anyway that the supposed speed advantage in going to 2^n points doesn't matter. And the result you want is badly affected.

Here is some code that 1) gets rid of nextpow2, 2) divides the fft by the number of points to scale the amplitude correctly, and 3) uses fftshift to put zero frequency in the middle, along with an appropriate frequency grid for that. You get a peak at +16.66 and a peak at -16.66, each of abs(height) = 10. This is appropriate for a sine wave of amplitude 20.

SampleInt = 0.00000004; % interval at which scope samples
ScopeFreq = 16.667; %[Hz] Reading from O-scope
timelength = 1:1:1499312;
Time = timelength.*SampleInt;
Max = 20;
Vv=Max*sin(2*pi*ScopeFreq*Time);
figure(2)
plot(Time,Vv)
xlabel('Time [Sec]')
grid on
% Now Check the FFT function using the reference sinewave
Lv=length(Vv)
Dtv = SampleInt;        % Get the sampling rate
Fsv = 1/Dtv; 
Xv=fft(Vv)/Lv;                    % divide by Lv for scaling
mxv = fftshift(abs(Xv));          % put zero frequency in the middle
fv =  (-Lv/2:Lv/2-1)*Fsv/Lv;      % corresponding frquency grid 
figure(3)
stem(fv,mxv)
title('FFT of Raw Test Data')
xlabel('Freq. [Hz]')
xlim([-100,100])

[Calvin Corey]

Thank you for the great explanation! This method ended up working and I correlated it to the FFT help website's method. Here's a test script using both your method and the FFT website's method:

%%%%%%%%%%%%%%%%%%%% Simple FFT Script %%%%%%%%%%%%%%%%%%%
clc;clear all;close all

SampleInt = 0.0000004;  % interval at which scope samples
ScopeFreq = 16.667;     %[Hz] Reading from O-scope
plotcycles = 7;         % how many cycles you want to plot in the time domain ( must be >= highest harmonic input)  
Max = 20;               % Peak of the fundamental waveform

timelength = 1:1:((plotcycles/ScopeFreq)/SampleInt);
Time = timelength.*SampleInt;

% waveform to be FFT'd
Vv=Max*sin(2*pi*ScopeFreq*Time) + Max*.05*sin(2*pi*ScopeFreq*3*Time)+ Max*.2*sin(2*pi*ScopeFreq*5*Time)+ Max*.2*sin(2*pi*ScopeFreq*7*Time);

%% Plot the time Signal
figure(1)
plot(Time,Vv)
xlabel('Time [Sec]')
grid on

%% Take the FFT of the time signal (single sided) using the FFT 
Dtv = SampleInt;        % Get the sampling rate
Fs = 1/Dtv;
L=length(Vv)
nfft = 2^nextpow2(length(Vv)) % Length of FFT
X=fft(Vv);
P2 = abs(X/L);
x_mag = P2(1:L/2+1);
f=Fs*(0:(L/2))/L;

%% STEM Plot the FFT
figure(2)
stem(f,x_mag)
title('FFT of Raw Test Data')
xlabel('Freq. [Hz]')
xlim([0,200])
grid on

%% Take the FFT using alternative method (double sided)

Lv=length(Vv)
Dtv = SampleInt;                  % Get the sampling rate
Fsv = 1/Dtv; 
Xv=fft(Vv)/Lv;                    % divide by Lv for scaling
mxv = fftshift(abs(Xv));          % put zero frequency in the middle
fv =  (-Lv/2:Lv/2-1)*Fsv/Lv;      % corresponding frquency grid 

%% STEM Plot the FFT (double sided)
figure(3)
stem(fv,mxv)
title('FFT of Raw Test Data')
xlabel('Freq. [Hz]')
xlim([-200,200])
grid on

Hopefully this helps someone else!

Cal