How do I perform fft for each column in a matrix?

Hello,

I am currently trying to perform fft to each column in two matrices (AD and C). In the below example, I have referenced columns 1 in each matrix, however when I attempt to plot my results, my peaks are always at 50Hz, not matter which column reference I enter. I have two methods of codes below, could you please tell me which is closer to the correct method, and what is going wrong in the examples? Thank you for your help!

Code 1

ADC3fft = fft(AD(:,1));
ADC3abs = abs(ADC3fft/L);
PADC3 = ADC3abs(1:L/2+1);
f = fs*(0:(L/2))/L;
PADC3(2:end-1)=2*PADC3(2:end-1);
 
ADC4fft = fft(AD(:,2));
ADC4abs = abs(ADC4fft/L);
PADC4 = ADC4abs(1:L/2+1);
f = fs*(0:(L/2))/L;
PADC4(2:end-1)=2*PADC4(2:end-1);
 
ADF3fft = fft(AD(:,3));
ADF3abs = abs(ADF3fft/L);
PADF3 = ADF3abs(1:L/2+1);
f = fs*(0:(L/2))/L;
PADF3(2:end-1)=2*PADF3(2:end-1)
 
figure(3)
subplot(2,1,1);
plot(f,PADC3)
 
subplot(2,1,2);
plot(f,PADC4)

Code 2
f = (0:length(AD)-1)*100/length(AD);
 
ADC3fft = fftshift(fft(AD(:,1)));
mADC3 = abs(ADC3fft);
 
CC3fft = fftshift(fft(C(:,1)));
mCC3 = abs(CC3fft);
 
figure(3)
subplot(2,1,1)
plot(f,mADC3)
 
subplot(2,1,2)
plot(f,mCC3) 

[dpb]

Description

example

Y = fft(X) computes the discrete Fourier transform (DFT) of X using a fast Fourier transform (FFT) algorithm.

• If X is a vector, then fft(X) returns the Fourier transform of the vector.
• If X is a matrix, then fft(X) treats the columns of X as vectors and returns the Fourier transform of each column.

Just take FFT() of the 2D array and use the results by column.

If the sampling rate is the same for each column, then all the rest is also the same for each column; the peaks will be where they will be.

[William Taylor]

Thanks for your help.

Would that therefore mean that 'code 1' above would be the best approach?

And also, could you please ellaborate on your last point above? Not sure I followed correctly.

Appreciate the advice

[dpb]

No. The better approach would be something like

YAD=fft(AD);                        % FFT by column of AD array
P2AD=abs(YAD/L);                    % 2-sided PSD
P1AD=P2AD(1:L/2+1,:);               % 1-sided
P1AD(2:end-1,:)=2*P1AD(2:end-1,:);  % scale, don't double DC and Fmax components
f = fs*(0:(L/2))/L;                 % compute frequency assuming all columns at same sample rate

Then, P1AD contains all the one-sided PSDs for the AD array assuming each column is an independent channel.

You can iterate over the number of channels to plot() each independently or use waterfall to compare on one axis or display however you choose.

Creating multiple variables with successive numeric suffixes is a sure sign one is not using MATLAB most effectively.

As far as what the frequency content of the signal is, the above will be correct for all channels/columns IFF all channels were sampled at the same frequency. If not, you need to use the proper sample rate by channel/column.

We have no way to know what those parameters were.

If all signals show a 50 Hz peak and if you're somewhere that that is the common electric power line frequency, that might be indication your signals are contaminated with a ground fault and what you're seeing is just the power line frequency. In the US that would be 60 Hz, not 50.

[William Taylor]

Hi dpb,

Thanks for the above answer, I have implemented the code and it is definitely making more sense to me now.

I have one follow up questions that maybe you could help with. P1AD is a '1 x n double' so how would this be iterated over to plot multiple graphs, is it is just 1-dimensional? All channels are sampled at the same frequency so f = fs*(0:(L/2))/L; will indeed fit all iterations.

I am in the UK so the point about mains electricity may make sense.

Thanks again

[dpb]

".... P1AD is a '1 x n double' ..."

No. P2 and P1 are both 2D arrays of M columns for however many columns are in the array AD

Note carefully those assignments to P1; there's a ",:" in the addressing expression to pull all columns over the subset of rows.

There was a brief time before I came back and fixed it that the post didn't have it there because I cut 'n pasted from the original code and didn't make the edit initially until I remembered later hadn't done.

So, if you picked up the code in that time interval, you may be missing that--if you do only have one column then either that particular input array had ony one column or the column addressing portion of the subscript is missing.

To iterate over columns, just something like

for i=1:size(P1AD,2)
  figure
  plot(f,P1AD(:,i))
  ...
end

is the simplest possible to write...and will create the crude plot for each as a new figure. Salt to suit...

[dpb]

It's always a potential problem that instrumentation may be picking up contamination from the power supply -- the 50 Hz is a telltale sign of that possibility. Need to check all grounds, etc., etc., ...

Answer 1

[Mathieu NOE]

hi

may I Jump in the conversation ?

there are so many questions in this forum dealing with fft , and many times it is not always clear for everybody how to get the best result according to what kind of signals are being analysed (random, harmonics, stationnary or not ...).

I like to propose this little demo based on pwelch and spegram to give another point of view and also let the people play with fft size, overlap , windows and see the effect on frequency resolution, amplitude resolution etc...

@ William maybe this can help you also in your work - enjoy it

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dummy signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs = 5000;
samples = 255001;
dt = 1/Fs;

t = (0:dt:(samples-1)*dt);
omega = 2*pi*(25+20*t);

sensordata = randn(size(t)) + 5*sin(omega.*t);  % signal = linear chirp + noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% FFT parameters
samples = length(sensordata);
NFFT = 512;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% option 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sensor_spectrum, freq] = pwelch(sensordata,w,NOVERLAP,NFFT,Fs);

figure(1),plot(freq,20*log10(sensor_spectrum));
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude (dB)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% option 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

saturation_dB = 80;  % dB range scale (means , the lowest displayed level is 80 dB below the max level)
% fmin = 1;
% fmax = Fs/2;

[sg,fsg,tsg] = specgram(sensordata,NFFT,Fs,w,NOVERLAP);  

% FFT normalisation and conversion amplitude from linear to dB peak
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% saturation to given dB range
mini_dB = round(max(max(sg_dBpeak))) - saturation_dB;
sg_dBpeak(sg_dBpeak<mini_dB) = mini_dB;

% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);axis('xy');colorbar('vert');
title(['Spectrogram  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

[William Taylor]

Hi there, many thanks for this. I will save this example to look at in the future!

Not sure if it can be directly applied to my current signal, as it is not generated as a sinusoid in the same method as above (it is infact a series of EEG readings)

Thanks again for the suggestion and I'll be sure to use this for some practice

[Mathieu NOE]

You re welcome

FYI , have you looked at : https://www.mathworks.com/matlabcentral/fileexchange/56415-eeglabdevelopers-eeglab

Answer 2

[Steven Lord]

If you want to make it explicit that you're operating down the columns, from the documentation page: "Y = fft(X,n,dim) returns the Fourier transform along the dimension dim. For example, if X is a matrix, then fft(X,n,2) returns the n-point Fourier transform of each row."