Fast Fourier Transform on multiple .csv files in one plot

Hello, I have 26 CSV files with 1001 time domain data points in each.I'm trying to perform an FFT on the data to view the frequency spectrum.

I have been able create an FFT of each file individually, but was hoping to do an FFT on all of the files in one go so that they all show up on the same plot.

I don't know how I would go about doing this, so I was hoping someone might be able to help me with this problem.

Best Regards

Robert

Answer 1

[Star Strider]

I would read them all in, each as its own cell array element (one for each file), then, since they are all the same lengths, concatenate them (most likely horizontally since I suspect they’re column vectors), using cell2mat or simple cell addressing to convert them from cells to a double matrix, then save all of them as the double matrix in a .mat file. For helpful details, see FAQ: How can I process a sequence of files?

Then load the .mat file and do the fft. It operates on columns by default. The code between the top two plot figures in this documentation for fft is everyting you need to know about analysing and plotting them.

[Robert Evans]

Thank you for your reply. I've been able to save all my data as a numeric matrix, 'data', and am attempting to perform the fft in one go using the code:

   N = length(data);

   freq = fft(fft(data))/N;

   plot(abs(freq));

However this provides a plot of waveforms rather than the spikes at specific frequencies that I would have expected i.e. I was expecting all the different FFT representations of my columns to line up in a row.

[Star Strider]

My pleasure.

It would, because taking the Fourier transform of Fourier transformed data is essentially the same as taking the inverse Fourier transform of it (but not actually, since the difference is in the sign of the argument in the complex exponential).

The code for taking the fft of your data is in the fft documentation that I provided here: fft. I’m copying it here, although you may need to alter it to fit your data:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 
y = x + 2*randn(size(t));     % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

The first block of code creates the signal. (You don’t need that section since you already have your signals. I included it because it explains some of the variables used in the rest of the code.) The second block of code does the actual fast-Fourier transform. The third block of code plots it.

You should be able to take the fft of your data in one go using an appropriately tweaked (for your data) version of the second and third blocks of the code here.

In my code, I create an additional variable, ‘Fn’, to designate the Nyquist frequency (the highest resolvable frequency in a sampled signal):

Fn = Fs/2;

Calculating it separately makes life easier for me. It’s just a personal preference. (You see it here included in the ‘f’ assignment.)

[Robert Evans]

Thanks again for your help. The plot shows the data on the x asis as an exponential i.e. it shows 2.5*10^8 instead of 250,000,000 (which would allow me to use hz as a label) or 250 ( which would alow me to use MHz as a label).

Do you know how I manipulate how the data on the axis is expressed?

[Star Strider]

I hope someday we will have the ability to format the tick labels without having to go through all the conniptions, but we’re not there yet. The only way to do that at present is to define the tick labels as strings.

Two examples:

x = linspace(0, 2.5E+8, 10);
y = randi(10, 1, 10);

figure(1)
plot(x, y)
grid
xlabel('Frequency (Hz)')

figure(2)
plot(x, y)
grid
xt = get(gca, 'XTick');
xts = regexp(sprintf('%5.1f\n', xt*1E-6), '\n', 'split');
set(gca, 'XTick',xt, 'XTickLabel',xts(1:end-1))
xlabel('Frequency (MHz)')

I included one digit to the right of the decimal in the '%5.1f' format string in the sprintf call. The field width of 5 allows three positions to the left of the decimal, the decimal point, and one digit to the right of the decimal. (You would also need to change it to 6 to accomodate a sign if there were negative values.) Change the format if you need to.

[Robert Evans]

Great, thanks very much for your help.

[Star Strider]

My pleasure.