MATLAB Answers

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

35 views (last 30 days)
William Taylor
William Taylor on 27 Oct 2020
Edited: dpb on 28 Oct 2020 at 15:39
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)

  6 Comments

Show 3 older comments
William Taylor
William Taylor on 28 Oct 2020 at 12:07
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
dpb on 28 Oct 2020 at 15:31
".... 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
dpb on 28 Oct 2020 at 15:34
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., ...

Sign in to comment.

Answers (2)

Mathieu NOE
Mathieu NOE on 27 Oct 2020 at 17:52
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)');

  2 Comments

William Taylor
William Taylor on 28 Oct 2020 at 12:09
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

Sign in to comment.


Steven Lord
Steven Lord on 27 Oct 2020 at 19:31
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."

  0 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by