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

5 次查看（过去 30 天）
Jake
Jake 2025-1-9
编辑： dpb 2025-1-12
Ran in:
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!
  5 个评论
显示 3更早的评论
dpb
dpb 2025-1-10
"... 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
Walter Roberson 2025-1-10
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

dpb
dpb 2025-1-11
编辑：dpb 2025-1-12
Ran in:
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')

更多回答（1 个）

Walter Roberson
Walter Roberson 2025-1-9
[~, 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.

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by