I'm trying to perform FFT on the 2nd data column of a .csv file. The .csv file is quite large and the zipped version can be found here: https://www.dropbox.com/s/87l5vtxdqfbtt6x/Tridral_1m_1ms.7z?dl=0
The signal is 10 milliseconds long, and is sampled 1 million times over that time period. I've tried following the FFT example provided by the help file, but the frequency domain signal does not come out right. The FFT output gives a strong DC signal (that could be right), and a signal at 50 KHz. Based on the time domain signal, i should expect approximately a 1 KHz fundamental with recurrent harmonics at 2, 3, 4 KHz..., but it comes out empty... am i doing something wrong in my code?

 Accepted Answer

dpb
dpb on 16 Feb 2019
Follow the recipe provided precisely...you've doubled the DC and Fmax components so that the DC bias is shown as twice what it would really be--and to find the frequency content in the signal you should either remove the mean first or, alternatively, just zero-out the DC component of the PSD. The former will be better numerically in preventing loss of precision in low-energy bins in the FFT, so that would be my recommendation.
...
y=raw_data(:,2)-mean(raw_data(:,2)); % eliminate DC bias
L=length(y);
Y = fft(y); % with 1E6 points, don't need interpolating
P2 = abs(Y/L); % two-sided spectrum
P1 = P2(1:L/2+1); % one-sided
P1(2:end-1) = 2*P1(2:end-1); % normalize; except not DC and Fmax that are only one point
f = Fs*(0:(L/2))/L;
plot(f,P1) % should show you actual frequency content
Undoubtedly, using semilogy() to amplify the dynamic range and zooming in on the lower frequency range will help to be able to see what structure there is in the signal.
You'll need to ensure the input signal is actually quality data and not contaminated somehow, too, of course...

4 Comments

Thanks for the tip! However, the output still seems odd. I've implemented the FFT part according to your inputs and the FFT outputs a fundamental at 10 Hz with harmonics at every 10 Hz (refer to '190217 FFT output_zoomed.png'. For a 1 ms periodic waveform, I'd expect a 1 KHz fundamental instead...is it something to do with my signal length or something?
After the FFT, i also find that it plots up to 500 KHz spectrum. Why is that so, when my signal is only about 1 KHz?
%setup signal for fft
Fs = 1000000; % sampling frequency in Hz
T = 1/Fs; % sampling period
L = 1000000; % signal length in ms
t = (0:L-1)*T; % time vector
% Read .csv, separate data from headers.
filename = 'Tridral_1m_1ms.csv'; % define filename
time_max = 10; % define horizontal time scale max(ms).
delimiterIn = ','; % define data delimiter
headerlinesIn = 21; % define header lines in .csv file.
raw = importdata(filename, delimiterIn, headerlinesIn); %import data
y = raw.data(:,2) - mean(raw.data(:,2)); % eliminate DC bias
% Plot imported data in time domain
datacount = size(raw.data,1); % obtain number of data points
timelength = [0:(time_max/datacount): (time_max - time_max/datacount)]; % create time X-axis
figure; hold on; grid on;
title('Signal in Time'); xlabel('Time (ms)'); ylabel('Voltage (V)');% setup figure
plotwave = plot(timelength,y); % plot time signal
L = length(y);
Y = fft(y); % with 1E6 points, don't need interpolating
P2 = abs(Y/L); % two-sided spectrum
P1 = P2(1:L/2+1); % one-sided
P1(2:end-1) = 2*P1(2:end-1); % normalize; except not DC and Fmax that are only one point
f = Fs*(0:(L/2))/L;
figure;
plot(f,P1);
dpb
dpb on 17 Feb 2019
Edited: dpb on 17 Feb 2019
"After the FFT, i also find that it plots up to 500 KHz spectrum."
Because you said you sampled at 1 MHz so Nyquist frequency is 1/2 Fmax --> 500 kHz
Also, from sampling relationships,
T = N dt and you said your T was 10 msec and had 1E6 samples.
So, the sampling RATE isn't 1 MHz, it's
1E6 samples/10 msec -->
1E6/(10*1E-3) --> 1E7 --> 10 MHz
So, you lied to the FFT in computing the frequency bins by having said you only sampled at a tenth of the speed you did--presuming the above numbers are correct, of course.
Set
Fs=10E6
instead and see if doesn't then make sense. Of course, the frequency resolution is also
df=1/T --> 1/(10*10E-3) = 10 Hz.
and Nyquist will be 5 MHz instead of 500 kHz.
Ray Fang
Ray Fang on 19 Feb 2019
Thanks. It makes sense now that the sampling rate is 1E6 samples/10ms = 100 MHz. I'm getting the fundamental at 1,KHz and harmonics at 2, 3..., KHz, which looks correct.
dpb
dpb on 19 Feb 2019
Yeah, I whiffed on that first go 'round; just took the Fs at face value without really looking at the rest that carefully...
Illustrates that the FFT doesn't care a whit about what the actual sample rate is; it just computes the numbers by bin; it's totally up to the user to set the scaling of the frequency axis to match the actual data collection process.

Sign in to comment.

More Answers (1)

Sofi Fatoni putra
Sofi Fatoni putra on 2 Oct 2019
Edited: Sofi Fatoni putra on 2 Oct 2019

0 votes

if I can ask for your data once again because the download in your dropbox can't be downloaded

Tags

Asked:

on 16 Feb 2019

Edited:

on 2 Oct 2019

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!