estimating oscillation frequency of a signal

I would like to find periodicity and the oscillation frequency in a signal. I used this simple MATLAB code which works find for some signals, but for a harmonic polluted signal like the one given below it generates wrong results. Is there any way to estimate the frequency of oscillation for a signal like the one given here. I appreciate any MATLAB code snippets.

[pks,pktimes] = findpeaks(x, (1:N)*DT);
Period = mean(diff(pktimes))

Answer 1

[Hassaan]

To estimate the frequency of a signal that is not perfectly harmonic (polluted with noise or other frequencies), you can use several approaches in MATLAB. One common method is to use the Fast Fourier Transform (FFT) to identify the dominant frequencies in the signal. Another method is to use signal processing functions like findpeaks to determine the period of the signal, from which the frequency can be calculated.

Generating Signal Data:

% MATLAB Code to Generate a Not Perfectly Harmonic Signal

% Sampling parameters
Fs = 1000;  % Sampling frequency
t = 0:1/Fs:1-1/Fs;  % Time vector

% Generate a dummy signal by combining different sine waves
f1 = 5; % Frequency of the first component (Hz)
f2 = 20; % Frequency of the second component (Hz)
signal = sin(2 * pi * f1 * t) + 0.5 * sin(2 * pi * 3 * f1 * t) + 0.25 * sin(2 * pi * 5 * f2 * t);
% signal = sin(2*pi*f1*t) + sin(2*pi*f2*t); % Signal with two frequency components

% Add noise to make the signal not perfectly harmonic
noise = 0.01 * randn(size(t));
signal_noisy = signal + noise;

% Plot the noisy signal
figure;
plot(t, signal_noisy);
title('Not Perfectly Harmonic Signal');
xlabel('Time (seconds)');
ylabel('Amplitude');
grid on;

Using FFT to Find Dominant Frequency:

% Assume 'signal_noisy' is your noisy signal and 'Fs' is the sampling frequency
N = length(signal_noisy); % Number of points in the signal
Y = fft(signal_noisy);    % Compute the FFT

% Two-sided spectrum
P2 = abs(Y/N);

% Single-sided spectrum
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);

% Frequency axis
f = Fs*(0:(N/2))/N;

% Find the peak in the spectrum
[~, idx] = max(P1);
dominantFrequency = f(idx);

% Plot the single-sided amplitude spectrum
figure;
plot(f, P1);
title('Single-Sided Amplitude Spectrum of the Signal');
xlabel('Frequency (Hz)');
ylabel('|P1(f)|');
grid on;

% Display the dominant frequency
disp(['Dominant frequency: ', num2str(dominantFrequency), ' Hz'])

Dominant frequency: 5 Hz

Note

The FFT method is more robust against noise, but it assumes that the signal is periodic within the window being analyzed.

-----------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

[Pooya Taheri]

@Muhammad Hassaan Shah

Thanks for your answers. The code you provided works. Here is the fundamental frequency plot of the data that I have which is quite accurate (0 most of the time, but jumps to 360, 60, 180 Hz at different times). How can I get rid of the glitches to have a smoother response?

Answer 2

[Sam Chak]

Hi @Pooya Taheri

The fft() command, or Fast Fourier Transform, can be used to estimate the oscillation frequency of a signal. When you apply the FFT to a time-domain sinusoidal signal, it transforms the signal from the time domain to the frequency domain, providing information about the various frequency components present in the signal. The resulting spectrum shows the magnitude and phase of each frequency component.

https://www.mathworks.com/discovery/fft.html

You can find various examples and code snippets in this documentation:

https://www.mathworks.com/help/matlab/ref/fft.html