How can I filter a FFt signal to only show the main frequency?

I have the following code (a sample data file is attached).

load("sample_signal.mat");
whos

  Name         Size             Bytes  Class     Attributes

  ans          1x41                82  char                
  signals      1x1             185952  struct              

figure
subplot(2,1,1)
hold on; grid minor
plot(signals.signal_01.signal)
plot(signals.signal_02.signal)
plot(signals.signal_03.time, signals.signal_03.theta, 'k')
title("Signal");

subplot(2,1,2)
hold on; grid minor
fs = 1;
[freq, amplitude, phase] = generic_fft(signals.signal_01.signal, fs);
plot(freq, amplitude)
[freq, amplitude, phase] = generic_fft(signals.signal_02.signal, fs);
plot(freq, amplitude)
fs = 10;
[freq, amplitude, phase] = generic_fft(signals.signal_03.theta, fs);
plot(freq, amplitude, 'k')
xlim([0 0.5])
xlabel("$\omega$ [rad/s]")
ylabel("fft")
title("FFT");

% fft function

function [freq, amplitude, phase] = generic_fft(signal, fs)
    % Remove NaNs
    signal = signal(~isnan(signal));
    % Length of the signal
    n = length(signal);
    % Perform FFT
    fft_values = fft(signal);
    % Compute single-sided amplitude spectrum
    amplitude = abs(fft_values / n); % Normalize by length
    amplitude = amplitude(1:floor(n / 2) + 1); % Keep positive frequencies
    amplitude(2:end-1) = 2 * amplitude(2:end-1);
    % Compute single-sided phase spectrum
    phase = angle(fft_values(1:floor(n / 2) + 1));
    % Frequency vector
    freq = (0:(n / 2)) * (fs / n);
end

Now, the signal plotted in black shows a clear, dominant frequency in the FFt, but the other two signals show additional "noise" (see, after 0.1 rad/s). Is it possible to maybe filter these out or only show the dominant frequency?

EDIT : What I wanted is not to identify one frequency point, but rather something similar to the black line in the fft plot. I.e., the signals used would need to be smoothened maybe?

The code is rough and the data structure is also not ideal (I saved the data mid-simulation for this question :) ). Also, any tips to improve the fft function is also very appreciated. TIA!

[Walter Roberson]

In an FFT plot, showing only the dominent frequency would be a pulse plot -- everything zero except for a single point.

[dpb]

Averaging will reduce any random components leaving the deterministic signals cleaner.

[Jake]

I agree with @Walter Roberson, and I should've worded my question better. My goal is not to get one frequency point, but to remove the noise from the signal so the fft would not show the noise as it is showing now. I'll edit the question here soon.

@dpb is right. I did try using a simple smooth command to smooth the signal before getting the fft, but it didn't get rid of the noise shown in fft. Maybe my fft is not right?

[dpb]

"... I did try using a simple smooth command to smooth the signal before getting the fft, ..."

Compute the PSD by averaging over sections of the signal to reduce stochastic noise in the result.

Subtract the mean or detrend to remove the DC component first...

[Walter Roberson]

It depends on the shape of your noise. For "salt and pepper" noise (gaussian noise), just zero some of the high frequency bins of the fft: gaussian noise is equivalent to adding a high frequency.

This will not have the effect of smoothing your fft output. Smoothing your fft output is equivalent to removing noise that is frequency-dependent.

Answer 1

[dpb]

Your signal is really quite short so it's hard to do too much about averaging other than using overlap, but the easy way of rounding the signal length to an even number and then reshaping and padding to the original length to keep the frequency resolution the same results for a case of half or averaging N=2 produces a noticeably smoother result already...this also did detrend to remove the DC and any slight trend that might be present.

load("sample_signal.mat");
figure
subplot(2,1,1)
hold on; grid minor
plot(signals.signal_01.signal)
plot(detrend(signals.signal_01.signal))
legend('Original','Detrended','location','south')
%plot(signals.signal_02.signal)
%plot(signals.signal_03.time, signals.signal_03.theta, 'k')
title("Signal");

subplot(2,1,2)
hold on; grid minor
fs = 1;
N=2;
[freq, amplitude, phase] = generic_fft(signals.signal_01.signal, fs);
plot(freq, amplitude)
[freq, amplitude, phase] = average_fft(signals.signal_01.signal, fs,N);
plot(freq, amplitude)
xlim([0 0.25])
xlabel("$\omega$ [rad/s]",'interpreter','latex')
ylabel("fft")
title("FFT");

function [freq, amplitude, phase] = generic_fft(signal,fs)
    % Remove NaNs
    signal = signal(~isnan(signal));
    n = length(signal);
    % Perform FFT
    fft_values = fft(signal);
    % Compute single-sided amplitude spectrum
    amplitude = abs(fft_values / n); % Normalize by length
    amplitude = amplitude(1:floor(n / 2) + 1); % Keep positive frequencies
    amplitude(2:end-1) = 2 * amplitude(2:end-1);
    % Compute single-sided phase spectrum
    phase = angle(fft_values(1:floor(n / 2) + 1));
    % Frequency vector
    freq = (0:fix(n/2))*fs/n;
end

function [freq, P1, phase] = average_fft(signal, fs,N)
    signal = signal(~isnan(signal));
    signal=detrend(signal);
    n = length(signal);
    M=N*fix(n/N);
    %L=N*round(ceil((n/N)),-2);
    signal=signal(1:M);
    signal=reshape(signal,[],N);
    signal=[signal;zeros(n-height(signal),N)];
    Y=mean(fft(signal),2);
    P2= N*abs(Y/n);
    P1=P2(1:floor(n/2)+1);
    P1(2:end-1) = 2*P1(2:end-1);
    % Compute single-sided phase spectrum
    phase = angle(Y(1:floor(n/2) + 1));
    % Frequency vector
    freq = (0:(n / 2)) * (fs / n);
end

If you have the Signal Processing TB, then this is easy with pwelch

figure, hold on
s=detrend(signals.signal_01.signal(isfinite(signals.signal_01.signal)));
[pxx,w] = pwelch(s);
plot(w/pi,pxx)
xlabel('\omega / \pi')
[pxx,w] = pwelch(s);
plot(w/pi,pxx)
N=2;
L=fix(numel(s)/N);
[pxx,w] = pwelch(s,L);
plot(w/pi,pxx)
N=4;
L=fix(numel(s)/N);
[pxx,w] = pwelch(s,L);
plot(w/pi,pxx)
legend('All','N/2','N/4')

Answer 2

[Walter Roberson]

[~, maxidx] = max(abs(amplitude));
main_signal = zeros(size(amplitude));
main_signal(maxidx) = amplitude(maxidx);

Do not be surprised if the peak frequency is maxidx == 1, corresponding to 0 Hz.