Frequency analysis using FFT

Dear colleagues, my question is How to analyse the frequencies and its corresponding magnitude within an array by using FFT function in MatLab?

I am certain that the values in the array have the following information:

• Already sampled at 1kHz
• Contains 30001 samples
• The sampled waveform is superimposed with multiple frequencies (e.g. 15Hz, 34Hz, etc)

I would like to analyse the magnitudes of each frequency band that is contained in the samples.

Thank you very much for any help!

Thomas

Answer 1

[Star Strider]

The fft documentation has everything you need in the code between the first two plot images. Your Nyquist frequency is 500 Hz.

To get the peak values and frequencies, use the Signal Processing Toolbox findpeaks function. It has a number of name-value pair arguments that make it extremely flexible.

[Thomas Leiw]

Thanks for your reply Star Strider.

I quickly did the following code with the attached CSV file:

filename = 'FFT_Case_1.4_A.csv';
signal = csvread(filename);

N = length(signal); % defines the total number of signal array
fax_bins = [0:N-1]; % fix the x-axis with the starting point of 0
half_N = ceil(N/2); % get the half number of samples for single side plot
fs = 1000; % sampling frequency of signal
fnyquist = fs/2; % Nyquist frequency of signal

signal_mag = abs(fft(signal)); % magnitude of signals after transform
fax_Hz = fax_bins*fs/N; % convert bins to Hz using simple formula

[PKS,FRQS] = findpeaks(signal_mag(1:half_N),'MINPEAKDISTANCE',200)

It seems that I need to re-index the magnitudes against frequencies instead of the sample number of the values. Would you mind to share with me the command that I could use for the re-indexing? Thanks.

Thomas

[Star Strider]

If you give findpeaks a second location vector argument (before the name-value pair arguments in the argument list), you can have the locations specified in terms of that vector rather than index values. See Find Peaks and Their Locations for details.

You apparently don’t specify a frequency vector in your code. If you define one, you can use the index locations as subscripts to it. You will need to define a frequency vector to use it either with subscripts or adding it to the findpeaks argument list.

Example:

d = xlsread('FFT_Case_1.4_A.csv');

Fs = 1E+3;
Fn = Fs/2;
Fd = fft(d)/length(d);                          
Fv = linspace(0, 1, fix(length(Fd)/2)+1)*Fn;    % Frequency Vector
Iv = 1:length(Fv);                              % Index Vector

[Pks,Frq] = findpeaks(abs(Fd(Iv)), Fv, 'MinPeakDistance', 20, 'MinPeakHeight',0.005);

figure(1)
plot(Fv, abs(Fd(Iv)))
hold on
plot(Frq, Pks, 'vr', 'MarkerFaceColor','m')
hold off
grid

[Thomas Leiw]

Sorry Star Strider, this might be a silly question ...

What is the reason for fft(d)/length(d)? Also, could you elaborate the line of

Fv = linspace(0, 1, fix(length(Fd)/2)+1)*Fn;    % Frequency Vector

Thanks.

Thomas

[Star Strider]

It’s necessary to normalise the fft so that the amplitudes are correct.

The ‘Fv’ assignment creates the frequency vector of half the length of the fft (to plot a one-sided fft) that goes from zero to ‘Fn’.

[Thomas Leiw]

Thanks for your explanation.

However, I am still not sure if I got the correct magnitude out of FFT. I am certain that the attached file contains the magnitude of almost 0.14 amps @15.8 Hz and 0.12 amps @16 Hz.

The following is the modified code:

filename = 'FFT_Case_1.4_A.csv';
signal = csvread(filename);

N = length(signal); % defines the total number of signal array
fax_bins = [0:N-1]; % fix the x-axis with the starting point of 0
half_N = ceil(N/2); % get the half number of samples for single side plot
fs = 1000; % sampling frequency of signal
fnyquist = fs/2; % Nyquist frequency of signal

signal_mag = abs(fft(signal)/(N/2)); % magnitude of signals after transform
fax_Hz = fax_bins*fs/N; % convert bins to Hz using simple formula

[PKS,FRQS] = findpeaks(signal_mag(1:half_N),'MINPEAKHEIGHT',0.1)

plot(fax_Hz(1:half_N), signal_mag(1:half_N))
title('Single-side Magnitude Spectrum')
xlabel('Frequency (Hz)')
ylabel('Magnitude (Amps)');

I divided the magnitude with the half of the signal length since only single side of the signal is at concern. Is there anything wrong in the code? Thanks a lot for your help.

Thomas

[Star Strider]

My pleasure.

I have a difficult time following your code. My code is a very slightly modified version of that in the fft documentation (I didn’t use nextpwr2, but it’s there simply to make the fft calculation more efficient), so I believe my code is correct.

I divided the entire signal by the entire length of the signal, then plotted only half of it to make it easier to interpret.

Your stated amplitude values may be correct, but your signal contains broadband noise that appears to have the characteristics of ‘1/f’ noise, and likely has energy at your frequencies of interest. That would result in a fft amplitude at those frequencies that is higher than what you may have measured. The fft has no way of distinguishing signal from noise at a particular frequency.

[Thomas Leiw]

The values that were mentioned in the previous post were measured magnitudes at those frequency bands using another software even if noise has contributed to it. I am curious why the fft function in MatLab returned different results. I drafted another code to demonstrate the magnitude of a superimposed signal as the following:

f1=14; % 1st frequency of 14 Hz
f2=40; % 2nd frequency of 40 Hz
fs=1000; % sampled frequency
L=1000; % length of the signal
T=1/fs; % sampled period
t=(0:L-1)*T; % total last time of the signal
signal = 0.7*sin(2*pi*f1*t)+sin(2*pi*f2*t); % superimposed sine waveforms

N = length(signal); % defines the total number of signal array
fax_bins = [0:N/2-1]; % fix the x-axis with the starting point of 0
half_N = ceil(N/2); % get the half number of samples for single side plot
fnyquist = fs/2; % Nyquist frequency of signal

signal_FFT = abs(fft(signal)/N);
% magnitude of single-sided signal after transform
signal_mag = 2*signal_FFT(1:N/2+1); 
fax_Hz = fax_bins*fs/N; % convert bins to Hz using simple formula

[PKS,FRQS] = findpeaks(signal_mag(1:half_N),'MINPEAKHEIGHT',0.69)

plot(fax_Hz, signal_mag(1:half_N))
title('Single-side Magnitude Spectrum')
xlabel('Frequency (Hz)')
ylabel('Magnitude (Amps)');

Since fft has the mirrored image at both sides of nyquist frequency, the magnitude needs to be multiplied by 2 in order to get the correct magnitude of 0.7 and 1 at the defined frequencies.

I think this is also implemented on the fft documentation page: http://au.mathworks.com/help/matlab/ref/fft.html

where:

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Do you think this is correct?

Thomas

[Star Strider]

It is correct, and entirely consistent with the previous version of the documentation on the fft function. (The documentation for fft changed in R2015b. The previous versions seemed to me to be more straightforward.) The single-sided fft in the documentation is correct, since it seems to have been taken from the R2015b fft documentation.

For your data, the amplitudes in my code at 16.3 Hz even with slightly revised code (corresponding to the R2015a fft documentation), are about 0.1 of those stated.

The only significant change in my code is:

[Pks,Frq] = findpeaks(2*abs(Fd(Iv)), Fv, 'MinPeakDistance', 0.25, 'MinPeakHeight',0.01);