By "middle" frequency, I'm assuming you want to filter out the signal so that only the component at 159 Hz is present. This means that the output result should only contain one sinusoid at 159 Hz. First off, I'd like to mention that the Fourier Transform shows a frequency decomposition of your signal. You can think of a signal as a summation of many sinusoids at different frequencies and they're possibly out of phase with different phase shifts. For each sinusoid at a frequency, there is a magnitude component and a phase component. The magnitude tells you how large the sinusoid is at this frequency, and the phase tells you how much of a delay that sinusoid is experiencing.
Now, when computing the FFT, this performs the Cooley-Tukey algorithm, which is a very efficient way of computing the Fourier Transform. It's computed in such a way where the first half of the signal shows the spectrum from 0 <= f < fn and the second half of the signal shows the spectrum from -fn <= f < 0. fn is what is known as the Nyquist frequency, which is half of the sampling frequency. Note the exclusivity in the range of the first half of the spectrum compared to the second half. We don't include fn in the first half, but we include -fn in the second half.
This only applies for real signals, which is your case with the three sinusoidal summed together. In order to make this a bit more meaningful, you can use fftshift so that you are reorganizing the spectrum so the 0 Hz / DC frequency appears in the middle, which will allow you to plot the spectrum so that it's between -fn <= f < fn.
Also, the way the signal is being plotted is normalized, which means that you see the frequencies are actually between -1 <= f <= 1. To translate this to actual frequencies, you use the following relationship:
freq = i * fs / N
i is the bin number or the point on the FFT you want, which in your case is from 0 up to 500. Remember, the first half of the signal represents the frequency distribution from 0 up to fn and from 0 up to 500 in the signal is what you need. The negative frequencies can be found by just substituting i for -i. fs is the sampling frequency and N is the size of the FFT, which in your case is the total length of the signal, 1001. To generate the right frequencies that correspond to the right points, you can use linspace and generate N+1 points between -fn and fn to ensure the spacing is correct, but because we don't include fn at the positive end, we remove this from the end of the range.
Also, to plot the magnitude and phase of the signal, use abs and angle respectively.
Therefore, try plotting it this way, and also focus on the magnitude and phase too. Just plotting the spectra is undefined because in general, there are real and imaginary components.
%// Your code
fs = 1000;
t = (0 : 1/fs : 1)';
f = [ 86, 159, 392 ];
x = sum( cos(2*pi * t * f), 2 );
y = fft(x);
%// New code
%// Shift spectrum
ys = fftshift(y);
N = numel(x); %// Determine number of points
mag = abs(ys); %// Magnitude
pha = angle(ys); %// Phase
%// Generate frequencies
freq = linspace(-fs/2, fs/2, N+1);
freq(end) = [];
%// Draw stuff
figure;
subplot(2,1,1);
plot(freq, mag);
xlabel('Frequency (Hz)');
ylabel('Magnitude');
subplot(2,1,2);
plot(freq, pha);
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');
Here's what I get:
As you can see, there are three spikes which correspond to the three sinusoidal components. You want the middle one, and that's at 159 Hz. To create a bandpass filter, you want to filter out all components except for the ones at +/- 159 Hz. If you want to do this automatically, you'd find the bin locations that are the closest to +/- 159 Hz, then expand a neighbourhood that surrounds these two points and ensure they're untouched while zeroing out the rest of the components.
Because you have exact sinusoids, using a bandpass filter in this regard is perfectly acceptable. In general, you don't do this due to ringing and aliasing effects as the sharpness of the cutoff from the bandpass filter in this way introduces unwanted aliasing effects in the time-domain. See the Wikipedia article on aliasing for more details.
So, to figure out where we need to filter out, try using min and find the absolute difference between the generated frequencies above with 159 Hz - specifically find the location. Once you find these for both +159 Hz and -159 Hz, extend a neighbourhood around these points, make sure they're untouched while the rest of the points are set to 0 in the spectrum:
[~,min_pt_pos] = min(abs(freq - f(2))); %// Find location where +159 Hz is located
[~,min_pt_neg] = min(abs(freq + f(2))); %// Find location of where -159 Hz is
%// Neighbourhood size
ns = 100; %// Play with this parameter
%// Filtered signal
yfilt = zeros(1,numel(y));
%// Extract out the positive and negative frequencies centered at
%// 159 Hz
yfilt(min_pt_pos-ns/2 : min_pt_pos+ns/2) = ys(min_pt_pos-ns/2 : min_pt_pos+ns/2);
yfilt(min_pt_neg-ns/2 : min_pt_neg+ns/2) = ys(min_pt_neg-ns/2 : min_pt_neg+ns/2);
yfilt now contains the filtered signal that removes all components except for 159 Hz. If you want to show the magnitude and phase of this filtered signal, we can do that:
mag2 = abs(yfilt);
pha2 = phase(yfilt);
figure;
subplot(2,1,1);
plot(freq, mag2);
xlabel('Frequency (Hz)');
ylabel('Magnitude');
subplot(2,1,2);
plot(freq, pha2);
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');
Here's what we get:
As we expect, there is only one strong frequency that decomposes this signal, and that's at 159 Hz. Now to reconstruct this signal, you'll have to undo the centering we did, and you'll have to use ifftshift on this filtered result, then do the inverse via ifft. You may also get small residual imaginary components, so it's a good idea to use real on the output result.
out = real(ifft(ifftshift(yfilt)));
If we plot this, we get:
plot(t, out);
xlabel('Time (seconds)');
ylabel('Height');
As you can see, there is a sinusoid with one frequency and that's at 159 Hz. Don't mind the amplitudes though. This is simply due to an imprecision in how many points you have chosen to plot the signal and so some of the time points may not exactly coincide with the true peaks of the signal. Do bear in mind that if there was more than one sinusoid present, there will be a point where the peaks of all sinusoids will meet and you'll get a higher peak instead of a peak of 1 that a single sinusoid offers. Because the amplitude hovers around -1 to 1, you can be certain that only one sinusoid is present. This can be avoided if you were to choose a finer grain step size, and hence more points in the FFT.
Hope this is enough to get you started. Good luck!
