Choosing one peak out of two in FFT function

Hello,

I have a question about how to choose only one out the two peaks that are produced after applying FFT function.

I have multiple text files where each one has "y" values, each file is at a certain time step so the peak of the FFT function is different for every file. Each one produces two peaks the same way in the image, is there a way to take just one of the highest two peaks for each file and store them to either plot or post process.

I tried findpeaks, but the problem is since it is varying peaks, it is difficult to set a parameter such as 'MinPeakDistance' or other.

I use the following:

A = fft(y);

B = abs(A);

figure()

hold on

semilogy(B)

hold off

even if I use C = max(abs(fft(y))); it would not decide which max to choose

Note: the first image is for one file, the second image is for multiple files

Thank you!

Answer 1

[William Rose]

@Hussein Kokash,

The output of y=fft(x) is symmetric about the middle of the array y.* This is due to the mathematics of the discrete Fourier transform. So when identifying peaks, you can ignore those that occur above element N/2, where N is the length of x and of y.

*When the input sequence,x, is real numbers.

[Hussein Kokash]

Any way to save every peak in different array? say peak1 and peak2 for each file, save in seperate arrays.

something like, P1 (saving the 1st peak for all files), P2 (saving the 2nd peak for all files).

The only reason is to extract the corresponding two x-positions (frequency) for each file

Thanks again William!

[William Rose]

@Hussein Kokash,

Here is an example using a waveform whose Fourier transform is similar to the one you showed. I am using your notation for y and B.

N=101;                      %number of points in signal
dt=1;                       %sampling interval (s)
t=(0:N-1)*dt;               %time values
y=sawtooth(2*pi*10*t/N,.3); %triangle wave
B=abs(fft(y));              %discrete Fourier transform
subplot(211); plot(y)       %plot y=original waveform
title('y')
subplot(212); plot(B)       %plot B=abs(fft(y))
title('B')

Since you want to find the peaks, I assume you also want to find the frequencies of the peaks. If so, you need to know the frequencies that go with the fft. Given sampling interval dt, and signal duration N (in points), the frequencies associated with fft(y) are

f=(0:N-1)/(N*dt);   %frequencies (in Hz, if dt is in seconds)

With this information, you can plot B versus frequency.

figure;
subplot(211); plot(t,y); 
title('y'); xlabel('Time (s)')
subplot(212); plot(f,B); 
title('B'); xlabel('Frequency (Hz)'); grid on

Since B is symmetric about its middle, let's keep only the lower half of B (and f):

N2=floor(N/2)+1; %number of points in first half of B, plus 1 
B=B(1:N2); f=f(1:N2);

Now do findpeaks on B to find the peak frequencies.

[pks,freqs] = findpeaks(B,f);
[pk1,idx1]=max(pks);        %biggest peak and its index
f1=freqs(idx1);             %frequency of biggest peak
hold on                     %hold existing plot
plot(freqs,pks,'g+')        %plot + at every peak
plot(f1,pk1,'r*')           %plot * at biggest peak

Try it.

[Image Analyst]

Call fftshift before plotting and it will make more sense.

[William Rose]

@Hussein Kokash, I had not read your comment when I posted my comment. Let me think about your question...

The code I gave in my comment makes a vector freqs, which has the frequencies of all the peaks in the lower half of the fft. And it makes a scalar f1 with the frequency of the biggest peak. Due to the mathematics of the discrete Fourier transform, the frequencies of the upper half are the mirror image of the lower half frequencies. The middle frequency is called the Nyquist frequency. It is exactly half the sampling frequency:

where dt is the sampling interval, in units of time. Therefore the frequencies of the upper half peaks are

fs=1/dt;  %dt was specified earlier in my code
freqs_up=fs-freqs;
f1_up=fs-f1;

Good luck.

[William Rose]

@Hussein Kokash, Why do you want to save the frequencies of the peaks in the upper half of the FFT? The reason I ask is that the upper half frequencies do not contain any additional useful informaiton.

Analogy: If you know and have saved p=probability of heads, for a particular coin, there is no need to save q=1-p=probability of tails.

[Hussein Kokash]

Hello William,

I appreciate your responces, so I am investigating a flow behind an object, I have a line of points behind this object (in the spanwise direction).

At each time step, I have a file that contains the spanwise (z) velocity for all of the points that makes this line (zmin to zmax), this is a sample:

0 -2.9577131E-12 0.236238

0.002 2.2367207E-10 0.236238

0.004 3.4525567E-10 0.236238

0.006 3.1392004E-10 0.236238

0.008 1.7650122E-10 0.236238

0.01 -7.5099878E-12 0.236238

0.012 -1.8817284E-10 0.236238

0.014 -3.1705425E-10 0.236238

0.016 -3.3854163E-10 0.236238

0.018 -2.1011342E-10 0.236238

0.02 1.8830366E-11 0.236238

0.022 2.4242510E-10 0.236238

0.024 3.5704751E-10 0.236238

0.026 3.1862220E-10 0.236238

0.028 1.7535199E-10 0.236238

0.03 -1.3538867E-11 0.236238

0.032 -1.9900941E-10 0.236238

0.034 -3.3320666E-10 0.236238

0.036 -3.6012169E-10 0.236238

0.038 -2.3574904E-10 0.236238

0.04 -7.6009011E-12 0.236238

0.042 2.1944382E-10 0.236238

0.044 3.4179600E-10 0.236238

0.046 3.1526771E-10 0.236238

0.048 1.8628746E-10 0.236238

0.05 1.0866790E-11 0.236238

0.052 -1.6580260E-10 0.236238

0.054 -2.9918222E-10 0.236238

0.056 -3.3471170E-10 0.236238

0.058 -2.2668630E-10 0.236238

0.06 -1.6748083E-11 0.236238

0.062 1.9794763E-10 0.236238

0.064 3.1805352E-10 0.236238

0.066 2.9838742E-10 0.236238

0.068 1.8193196E-10 0.236238

0.07 2.0337331E-11 0.236238

0.072 -1.4520383E-10 0.236238

0.074 -2.7288414E-10 0.236238

0.076 -3.0918668E-10 0.236238

0.078 -2.0794240E-10 0.236238

0.08 -9.7146442E-12 0.236238

0.082 1.9173136E-10 0.236238

0.084 3.0145340E-10 0.236238

0.086 2.7771161E-10 0.236238

0.088 1.6262074E-10 0.236238

0.09 5.5778040E-12 0.236238

0.092 -1.5427218E-10 0.236238

0.094 -2.7635975E-10 0.236238

0.096 -3.0763103E-10 0.236238

0.098 -2.0193161E-10 0.236238

0.1 3.3841662E-13 0.236238

0.102 2.0144567E-10 0.236238

0.104 3.0403225E-10 0.236238

0.106 2.6902292E-10 0.236238

0.108 1.4431987E-10 0.236238

0.11 -1.6204361E-11 0.236238

0.112 -1.7232844E-10 0.236238

0.114 -2.8591619E-10 0.236238

0.116 -3.0836106E-10 0.236238

0.118 -1.9723919E-10 0.236238

0.12 7.6309077E-12 0.236238

0.122 2.1455048E-10 0.236238

0.124 3.2751216E-10 0.236238

0.126 3.0043147E-10 0.236238

0.128 1.7559393E-10 0.236238

0.13 4.9337891E-12 0.236238

0.132 -1.6878562E-10 0.236238

0.134 -3.0112440E-10 0.236238

0.136 -3.3568558E-10 0.236238

0.138 -2.2463601E-10 0.236238

0.14 -1.0262878E-11 0.236238

0.142 2.1156070E-10 0.236238

0.144 3.3978994E-10 0.236238

0.146 3.2244513E-10 0.236238

0.148 1.9825914E-10 0.236238

0.15 1.9627189E-11 0.236238

0.152 -1.6648793E-10 0.236238

0.154 -3.0945631E-10 0.236238

0.156 -3.4744837E-10 0.236238

0.158 -2.3123037E-10 0.236238

0.16 -9.0491139E-12 0.236238

0.162 2.1339540E-10 0.236238

0.164 3.3424095E-10 0.236238

0.166 3.0936684E-10 0.236238

0.168 1.8137765E-10 0.236238

0.17 3.1284328E-12 0.236238

0.172 -1.8034544E-10 0.236238

0.174 -3.2127409E-10 0.236238

0.176 -3.5975507E-10 0.236238

0.178 -2.4529187E-10 0.236238

0.18 -1.9049244E-11 0.236238

0.182 2.1396086E-10 0.236238

0.184 3.4508396E-10 0.236238

0.186 3.2556036E-10 0.236238

0.188 1.9843574E-10 0.236238

0.19 1.9107855E-11 0.236238

0.192 -1.6491127E-10 0.236238

0.194 -3.0476838E-10 0.236238

0.196 -3.4132776E-10 0.236238

0.198 -2.2626663E-10 0.236238

0.2 -2.9577131E-12 0.236238

1st column (z), 2nd column (coresspoding z-velocity), 3rd colum (corresponding time)

The idea is I have say many files like these, where time step is variant.

The reason that I want to know the frequency at of each peak for each file is to see if there is a change in wavelength in time.

What makes it challenging is there are two main frequencies and if I just used the regular max or findpeaks, it would not be quite accurate which one was captured as a peak/frequency.

Thanks again!

[William Rose]

@Hussein Kokash,

You're welcome.

You have an interesting problem and an interesting experiment.

The code I provided gives f1, the frequency with the most power in the signal. Obviously the associated wavelength is lambda1=c/f1. The frequencies of the peaks in the "top half" of the spectrum (the part above the Nyquist frequency) have no additional value for your analysis.

Good luck with your work!